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(57) ABSTRACT 

A system and method for detecting structural damage is 
provided that utilizes a general order perturbation method- 
ology involving multiple perturbation parameters. The per- 
turbation methodology is used iteratively in conjunction 
with an optimization method to identify the stiffness param- 
eters of structures using natural frequencies and/or mode 
shape information. The stiffness parameters are then used to 
determine the location and extent of damage in a structure. 
A novel stochastic model is developed to model the random 
impact series produced manually or to generate a random 
impact series in a random impact device. The random impact 
series method or the random impact device can be used to 
excite a structure and generate vibration information used to 
obtain the stiffness parameters of the structure. The method 
or the device can also just be used for modal testing 
purposes. The random impact device is a high energy, 
random, and high signal-to-noise ratio system. 
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SYSTEM AND METHOD FOR DETECTING 
STRUCTURAL DAMAGE 

[0001] This application churns priority to U.S. Provisional 
Patent Application Ser. No. 60/471,873 filed May 20, 2003 
and U.S. Provisional Patent Application Ser. No. 60/512,656 
tiled Oct. 10, 2003, which arc both incorporated herein by 
reference. 

BACKGROUND OF THE INVENTION 
[0002] 1. Field of the Invention 

[0003] The invention relates to a method and apparatus for 
detecting structural damage, and, more specifically, to a 
method and apparatus for detecting structural damage using 
changes in natural frequencies and/or mode shapes. 
[0004] 2. Background of the Related Art 
[0005] Damage in a structure can be defined as a reduction 
in the structure's load bearing capability, which may result 
from a deterioration of the structure's components and 
connections. All load bearing structures continuously accu- 
mulate structural damage, and early detection, assessment 
and monitoring of this structural damage and appropriate 
removal from service is the key to avoiding catastrophic 
failures, which may otherwise result in extensive property 
damage and cost. 

[0006] A number of conventional non-destructive test 
(NDT) methods are used to inspect load bearing structures. 
Visual inspection of structural members is often unquanti- 
liable ami unreliable, especially in instances where access to 
damaged areas may be impeded or damage may be con- 
cealed by paint, rust, or other coverings. Penetrant testing 
(PT) requires that an entire surface of the structure be 
covered with a dye solution, and then inspected. PT reveals 
only surface cracks and imperfections, and can require a 
large amount of potentially hazardous dye be applied and 
disposed of. Similarly, magnetic particle testing (MT) 
requires that an entire surface of the structure be treated, can 
be applied only to ferrous materials, and detects only rela- 
tively shallow cracks. Further, due to the current required to 
generate a strong enough magnetic field to detect cracks, MT 
is not practically applied to large structures. Likewise, eddy 
current testing (ET) uses changes in the flow of eddy 
currents to detect flaws, and only works on materials that are 
electrically conductive. Ultrasonic testing (UT) uses trans- 
mission of high frequency sound waves into a material to 
detect imperfections. Results generated by all of these 
methods can be skewed due to surface conditions, and 
cannot easily isolate damage at joints and boundaries of the 
structure. Unless a general vicinity of a damage location is 
known prior to inspection, none of these methods are easily 
or practically applied to large structures which are already in 
place and operating. On the other hand, resonant inspection 
methods are not capable of determining the extent or loca- 
tion of damage, and is used only on a component rather than 
a assembled structure. None of the above NDT methods are 
easily or practically applied to large structures requiring a 
high degree of structural integrity. 

[0007] Because of these shortfalls in existing NDT meth- 
ods when inspecting relatively large structures, structural 
damage detection using changes in vibration characteristics 
has received much attention in recent years. Vibration based 
health monitoring for rotating machinery is a relatively 



mature technology, using a non-model based approach to 
provide a qualitative comparison of current data to historical 
data. However, this type of vibration based damage testing 
does not work for most structures. Rather, vibration based 
damage detection for structures is model based, comparing 
test data to analytical data from finite element models to 
detect the location(s) and extent of damage. Vibration based 
damage detection methods fall into three basic categories. 
The first of these is direct methods such as optimal matrix 
updating algorithms, which identify damage location and 
extent in a single iteration. Because of the single iteration, 
these methods are not accurate in detecting a large level of 
damage. The second category is iterative methods. The 
methodology has only been for updating modeling, which 
determines modified structural parameters iteratively by 
minimizing differences between model and test data. The 
third category includes control-based eigenstructure assign- 
ment methods, which have the similar limitation to that of 
the direct methods indicated above and are not accurate in 
detecting a large level of damage. None of these current 
v ibration based methods have been incorporated into an 
iterative algorithm that can detect small to large levels of 
damage, and the vibration based approach for structures 
remains an immature technology area which is not readily 
available on a commercial basis. 

SUMMARY OF THE INVENTION 

[0008] An object of the invention is to solve at least the 
above problems and/or disadvantages and to provide at least 
the advantages described hereinafter. 

[0009] Another object of the invention is to provide a 
system and method for delecting structural damage based on 
changes in natural frequencies and/or mode shapes. 
[0010] An advantage of the system and method as embod- 
ied and broadly described herein is that it can be applied to 
a large operating structure with a large number (thousands or 
more) of degrees of freedom. 

[0011] Another advantage of the system and method as 
embodied and broadly described herein is that it can accu- 
rately detect the location(s) and extent of small to large 
levels of damage and is especially useful for detecting a 
large level of damage with severe mismatch between the 
eigenparameters of the damaged and undamaged structures. 
[0012] Another advantage of the system and method as 
embodied and broadly described herein is that it can work 
with a limited number of measured vibration modes. 

[0013] Another advantage of the system and method as 
embodied and broadly described herein is that it can use 
measurement at only a small number of locations compared 
to the degrees of freedom of the system. A modified eigen- 
vector expansion method is used to deal with the incomplete 
eigenvector measurement problem arising from experimen- 
tal measurement of a lesser number of degrees of freedom 
than that of the appropriate analytical model. 
[0014] Another advantage of the system and method as 
embodied and broadly described herein is that it can be 
applied to structures with slight nonlinearities such as open- 
ing and closing cracks. The random shaker test or the 
random impact series method can be used to average out 
slight nonlinearities and extract linearized natural frequen- 
cies and/or mode shapes of a structure. 
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[0015] Another advantage of the system and method as 
embodied and broadly described herein is that it can handle 
structures with closely spaced vibration modes, where mode 
switching can occur in the damage detection process. 

[0016] Another advantage of the system and method as 
embodied and broadly described herein is that it can handle 
different levels of measurement noise with estimation errors 
within the noise levels. 

[0017] Another advantage of the system and method as 
embodied and broadly described herein is that the damage 
detection method and the vibration testing methods such as 
the random impact series method enables damage detection 
and assessment to be automated, thus improving the reli- 
ability integrity of results. 

[0018] Another advantage of the system and method as 
embodied and broadly described herein is that damage 
detection and assessment may be automated in the field so 
that structural health can be monitored at central location 
and useful service life may be optimized. 

[0019] Another advantage of the system and method as 
embodied and broadly described herein is that the random 
impact series method enables the modal parameters such as 
natural frequencies and/or mode shapes to be measured for 
a large structure or a structure in the field when there are 
noise effects such as those arising from the wind or other 
ambient excitation. 

[0020] Additional advantages, objects, and features of the 
invention will be set forth in part in the description which 
follows and in part will become apparent to those having 
ordinary skill in the art upon examination of the following 
or may be learned from practice of the invention. The objects 
and advantages of the invention may be realized and attained 
as particularly pointed out in the appended claims. 

BRIEF DESCRIPTION OF THE DRAWINGS 

[0021] The invention will be described in detail with 
reference to the following drawings in which like reference 
numerals refer to like elements wherein: 

[0022] FIG. 1A shows a system 100 for detection of 
structural damage according to one embodiment of the 
invention; 

[0023] FIG. IB is a flowchart of an inverse algorithm for 
identifying stiffness parameters of a damaged structure from 
a select set of measured eigenparamelers. in accordance with 
an embodiment of the invention; 

[0024] FIG. 2 is a flowchart of a quasi-Newton method for 
finding an optimal solution to the system equations shown in 
the flowchart of FIG. IB, in accordance with an embodi- 
ment of the invention; 

[0025] FIG. 3 is a schematic view of a serial mass-spring 

system; 

[0026] FIGS. 4A-4C illustrate estimation errors in a series 
of iterations for a low order system with a small level of 

damage; 

[0027] FIGS. 5A-5C illustrate estimation errors in a series 
of iterations for a low order system with a large level of 
damage; 



[0028] FIGS. 6A-6C illustrate estimation errors in a series 
of iterations for a large order system with a large level of 
damage; 

[0029] FIG. 7 illustrates a finite element model of a 
fixed-fixed beam; 

[0030] FIGS. 8A-8C illustrate the estimated stiffness 
parameters with complete eignenvector measurements and 
different noise levels for a ten element beam with a large 
level of damage; 

[0031] FIGS. 9A-9B illustrate stiffness parameters with 
reduced eigenvector measurements and different noise lev- 
els for a ten element beam with a large level of damage; 

[0032] FIG. 10 illustrates a modular, four bay space 
frame; 

[0033] FIG. 11 illustrates the estimated dimensionless 
stiffnesses for the damaged frame as shown in FIG. 11; 

[0034] FIG. 12 illustrates a cantilever aluminum beam test 
specimen with uniform damage in approximately 5 elc- 

[0035] FIG. 13 illustrates the experimental results for the 
estimated bending stiffnesses of all the elements of a can- 
tilever aluminum beam test specimen with the actual dam- 
age between 10 cm and 15 cm from the cantilevered end 
using a 40-element finite element model; 

[0036] FIG. 14 illustrates the experimental results for the 
estimated bending stiffnesses of all the elements of a can- 
tilevered aluminum beam test specimen with the actual 
damage between 25 cm and 30 cm from the cantilevered end 
using a 40-element finite element model; 

[0037] FIG. 15 illustrates the experimental results for the 
estimated bending stiffnesses of all the elements of an 
undamaged cantilever aluminum beam test specimen using 
a 40-element finite element model; 

[0038] FIG. 16 illustrates a cantilever aluminum beam test 

specimen with a narrow cut; 

[0039] FIG. 17 illustrates the experimental results for the 
estimated bending stiffnesses of all the elements of the test 
specimen as shown in FIG. 16 using a 45-element finite 
element model; 

[0040] FIG. 18 illustrates the estimated bending stiff- 
nesses of all the elements of the cantilever aluminum with 
simulated damage between 9 cm and 15.75 cm from the 
cantilevered end using a 40-element finite element model; 

[0041] FIG. 19 illustrates the estimated bending stiff- 
nesses of all the elements of the cantilever aluminum beam 
with simulated damage between 29.25 cm and 36 cm from 
the cantilevered end using a 40-element finite element 
model; 

[0042] FIG. 20 illustrates the estimated bending stiff- 
nesses of all the elements of the cantilever aluminum beam 
with multiple simulated damage — one at the 3 rd element and 
the other at the 20 th element; 

[0043] FIG. 21 illustrates the estimated bending stiff- 
nesses of all the elements of the cantilever aluminum beam 
with simulated uniform damage from the 16 th to the 18 th 
element; 
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[0044] FIG. 22 illustrates a 50-foot lightning mast test 
specimen with an eccentric spike at the top; 

[0045] FIG. 23 illustrates the experimental results for the 
estimated bending stiffnesses of all the elements of the 
lightning mast as shown in FIG. 22; 
[0046] FIG. 24 illustrates the estimated bending stiff- 
nesses of all the elements of the lightning mast as shown in 
FIG. 22 with simulated damage between 0.76 m and 1.125 
m from the ground using a 40-element finite element model; 
[0047] FIG. 25 illustrates the estimated bending stiff- 
nesses of all the elements of the lightning mast as shown in 
FIG. 42 with simulated damage between 6.89 m and 7.35 m 
from the ground using a 40-element finite element model; 
[0048] FIGS. 26A-26C illustrate the estimated bending 
stiffnesses of all the elements of a cantilever aluminum beam 
using the first regularization method; 

[0049] FIGS. 27A-27C illustrate the estimated bending 
stiffnesses of all the elements of a cantilever aluminum beam 
using the second regularization method; 
[0050] FIG. 28 illustrates an example of a practical appli- 
cation of the damage detection method of FIG. 1A, in 
accordance with an embodiment of the invention; 

[0051] FIG. 29 illustrates the application of a series of 
random impacts to the practical application of the damage 
detection method of FIG. 1A, in accordance with an 
embodiment of the invention; 

[0052] FIG. 30 illustrates a system for applying the series 
of random impacts shown in FIG. 29 to the practical 
application of the damage detection method of FIG. 1A, in 
accordance with an embodiment of the invention; 

[0053] FIG. 31 is a block diagram of one preferred 
embodiment of the random impact device of FIG. 30; 

[0054] FIG. 32 illustrates a lightning mast specimen; 
[0055] FIGS. 33A-33B are graphical representations of 
coherence and FRF from single and multiple impact tests of 
the mast shown in FIG. 32; 

[0056] FIG. 34 illustrates a random series of pulses with 
the same deterministic shape and random amplitudes and 
arrival times; 

[0057] FIG. 35 illustrates an average normalized shape 
function of force pulses; 

[0058] FIGS. 36-38 are graphical representations of ana- 
lytical and numerical solutions; 

[0059] FIG. 39 shows the comparison of the probability 
density function (PDF) of the experimental arrival time with 
that from uniform distribution; 

[0060] FIG. 40 shows the comparison of the experimental 
and analytical PDFs for the number of arrived pulses; and 

[0061] FIG. 41 shows the comparison of the experimental 
and analytical PDFs for the pulse amplitudes. 

DETAILED DESCRIPTION OF PREFERRED 
EMBODIMENTS 

[0062] Commonly measured modal parameters, such as 
natural frequencies and mode shapes, are functions of physi- 



cal properties of a particular structure. Therefore, changes in 
these physical properties, such as reductions in stiffness 
resulting from the onset of cracks or a loosening of a 
connection, will cause detectable changes in these modal 
parameters. Thus, if the changes in these parameters are 
indicators of damage, vibration based damage detection may 
be, simplistically, reduced to a system identification prob- 
lem. However, a number of factors have made vibration 
based damage detection difficult to implement in practice in 
the past. 

[0063] The system and method for detecting structural 
damage as embodied and broadly described herein is moti- 
vated by the observed advantages of vibration based damage 
detection over currently available technologies. It is well 
understood that this system and method may be effectively 
applied to damage detection and assessment for substantially 
all types and configurations of structures, including, but not 
limited to, simple beams, hollow tubes, trusses, frames, and 
the like. However, simply for ease of discussion, the system 
and method will first be discussed with respect to three 
examples — a mass-spring model, a beam, and a space 
frame — for conceptualization purposes. The system and 
method will later he applied lightning masts in electric 
substations. 

[0064] FIG. 1A shows a system 100 for detecting struc- 
tural damage according to one embodiment of the invention. 
System 100 includes a stiffness parameter unit 103 for 
detecting stiffness of a structure in question. Stiffness param- 
eter unit 103 may include a spectrum analyzer or any other 
device capable of receiving vibration information and pro- 
viding natural frequency and/or mode shape data. One 
example might be a spectrum analyzer capable of perform- 
ing a fast fourier transform (FFT) and with modal analysis 
software for deriving mode shape data from vibration infor- 
mation received from sensor 110 if mode shape information 
is needed. A roving sensor technique or multiple sensors 
may need to be used if mode shape information needs to be 
used. If mode shape information is desired, a sensor should 
he provided at the tip of the hammer or device used to excite 
the structure. This sensor should be configured to send data 
to the spectrum analyzer in order to obtain the frequency 
response function, which the spectrum analyzer (modal 
analysis software) uses to obtain mode shape information. 
Stiffness parameter unit 103 is configured to receive vibra- 
tion information output from the sensor 110 via sensor 
coupler 113 and input 114. Sensor coupler 113 could be a 
wire, optical fiber or even a wireless connection between 
sensor 110 and stiffness parameter unit 103. 

[0065] Stiffness parameter unit 103 may include an itera- 
tive processing unit 115 capable of determining stiffness 
parameters using a first order perturbation approach and the 
generalized inverse method. The first order perturbation 
approach can include using natural frequencies and/or mode 
shapes of the structure according to a preferred embodiment 
of the invention. Iterative processing unit 115 may be 
capable of converging to a correct result for the output of 
stiffness parameters even when the system is underdeter- 
mined. Stiffness parameter unit 103 may further include an 
outer iterative processing unit 117 and an inner (nested) 
iterative processing unit 119 which may operate using a first 
or higher order perturbation approach and a gradient or 
quasi-Newton method to be discussed below. Stiffness 
parameter unit 103 outputs stiffness parameters to damage 
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information processor 123. Damage information processor 
includes a damage extent processor 125 and a damage 
location processor 127. Damage location processor 127 
outputs the location(s) of damage and damage extent pro- 
cessor 125 outputs the extent of damage. 

[0066] System 100 using the iterative processing units 115 
or 117 and 119 is capable of determining the stiffness 
parameters and ultimately the location(s) and extent of 
damage for structures with a largest dimension less than 50 



, 25 r 



, 10 r 



meters. Examples of this are provided herein. 

[0067] The system and method may include a multiple- 
parameter, general order perturbation method, in which the 
changes in the stiffness parameters are used as the pertur- 
bation parameters. By equating the coefficients of like -order 
terms involving the same perturbation parameters in the 
normalization relations of eigenvectors and the eigenvalue 
problem, the perturbation problem solutions of all orders 
may be derived, and the sensitivities of all eigenparameters 
may be obtained. The perturbation method may be used in 
an iterative manner with an optimization method to identify 
the stiffness parameters of structures. 

[0068] Methodology 

[0069] This method presented below can simultaneously 
identify all l he unknow n stillness parameters and is formu- 
lated as a damage detection problem. Since the effects of the 
changes in the inertial properties of a damaged structure are 
usually relatively small, only the changes in the stiffness 
properties due to structural damage are considered. 

[0070] Consider a N degree-of-freedom, linear, time-in- 
variant, self-adjoint system with distinct eigenvalues. The 
stiffness parameters of the undamaged structure are denoted 
by G h; (i=l, 2, . . . , m), where m is the number of the stiffness 
parameters. Structural damage is characterized by reduc- 
tions in the stiffness parameters. The estimated stiffness 
parameters of the damaged structure before each iteration 

arc denoted by Ci ; (i=l. 2 in), and its stiffness matrix. 

which depends linearly on G ; , is denoted by K=K(G), where 
G^G-p G 2 , . . . , G m ] T . Here the superscript T denotes matrix 
transpose. The eigenvalue problem of the structure with 
stiffness parameters G ; is 



K$ k =X k M<|) k 



(1) 



[0071] where M is the constant mass matrix, and A. k =?i k (G) 
and <|> k =<|> k (G) (k=l, 2, . . . , N) are the k-th th eigenvalue and 
mass-normalized eigenvector respectively. It is noted that 
X k =o> k 2 , where w k is the k-th natural frequency of the 
structure. The normalized eigenvectors of (1) satisfy the 
orthonormality relations 



(<|) k ) T M<|)°=8 k u(1> k ) T A'q) u =>. k & liu 



(2) 



[0072] where l^u^N and 8^ is the Kronecker delta. 
Before the first iteration, the initial stiffness parameters of 
the damaged structure are assumed to be G ; (o:) =a ; G h; (i=l, 2, 
. . . , m), where 0<a ; s=l, and the eigenvalue problem in (1) 
corresponds to that of the structure with stiffness parameters 
G ; IU) . If there is no prior knowledge of the integrity of the 
structure, one can start the iteration from the stiffness 
parameters of the undamaged structure and set ] ; =1. Let G d; 
(i=l, 2, . . . , m) denote the stiffness parameters of the 



Kjtt-kfM*f (3) 
[0073] where K^I^GJ is the stiffness matrix with G d = 
[G d i,Gd2, • • • , G d J T , and X d k =X\G d ) and (|> d k =<|> k (G d ) are 
the k-th eigenvalue and mass-normalized eigenvector 
respectively. The stiffness matrix K d is related to K through 
the Taylor expansion 



[0074] where 8G ; =G di -Gj (i=l, 2, . . . , m) are the changes 
in the stiffness parameters, and the higher-order derivatives 
of K with respect to G ; vanish because K is assumed to be 
a linear function of G ; . Based on the finite element model, 
the global stiffness matrix of a distributed structure satisfies 
(4) as its higher-order derivatives with respect to each 
element stiffness parameter vanish. 

[0075] Let the k-th eigenvalue and mass-normalized 
eigenvector of the damaged structure be related to X k and <|> k 
through 



[0076] where k , ^ (2 )ij k > • • • , and X (p);j k are the 
coefficients of the first, second, . . . , and p-th order 
perturbations for the eigenvalue, z (1); k , z (2);j k , . . . , and z (p)ij 
. . . t k are the coefficient vectors of the first, second, . . . , and 
p-th order perturbations for the eigenvector, and e x k and e^ k 
are the residuals of order p+1. Note that the numbers in the 
parentheses in the subscripts of the coefficients and coeffi- 
cient vectors indicate the orders of the terms. By the Taylor 
expansion, one has for any p = l, 



[0077] By (7), 
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[0078] are symmetric in the p indices, i,j, . . . , and t. The 
right-hand sides of (7) are the p-th order sensitivities of the 
eigenvalues and eigenvectors with respect to the stiffness 
parameters. 



[0079] Using the normalization relations of the eigenvec- 
tors, (j) k and (|) d k , and symmetry of the coefficient vectors in 
(6), as indicated earlier, one obtains 



1 = (<fo T Mfc < 8 ' 



l + g[(/) T M4,. + (4 1)i ) r M/] 
dG, + 

XX jf^ )TM ^ + [ ( 4/ M *W + $u/ M &»] + 2'.(^ )iJ ) T M^}6G i SGj + 



XXX !</ )7mz " 3,! " + 2 'l^ 1 *^ 2 '--' + (4A4w. + (4/^ 2W ] + 



[(4 ) ,) 7 'm4_ 1)j .,... s , + (4„,i'a'4 , k , » + ■■■ + (4d/ m 4-iw J + 

2!(p - 2)![(4/m4_ 2) ,..„ + (4 2 ,-,) 7 'M^ 2)J -... s , + ... + (s^/Ms^,.. r ] + ... + 
(p -2)!2![(s^ 2y ... yAft^. + (^ 2)J ,.. s ,) 7 M z { 2 , i + .- + (4_ 2>y ... //W^J + 

(p - wM P - l)Jl ...s,) T Mz k a)i + (Ap-m-J M Ai)j + ■■■ + (4-i)i,.yMzy + 
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[0080] where the superscript T denotes transpose. For any 
p-th (p = l) order term in the last expression of (8), the 
coefficient 



[0081] is defined by 



4..., = X 1 !X 2 !-X 0 !, 



-continued 

(42)i2) 7 M4 2)22 + (4 22 /m4 12 + 

(42) 22) 7 M42)12 + (42>22) r <W4 12 + 

3![(4)222) 7 'Mzfi.i + <4:/ M 4: + 

(43) 122) ? 'm4, )2 + (4„ 22 ) 7 'Mzf 1)2 ] + 

4! (4t ) i22 2 ) T ' w ^) lSG i lSG : = 
{4(/) T M4 4) , 222 + [u-f,,, 1^4™ + ^4ic I r «4i::] + 

4(42)12) 7 ''W422 + [(4)222) T M4l + 

3l4 122 ) 7 'A74 2 ]+4(4 4)1222 /M^j I 5G 1I 5G| 



[0082] where X 1; X 2 , . . . , and X, (liaSp) are the 
numbers of the first, second, . . . , and last distinct indices 
within indices, i, j, . . . , and t, with Xj+X 2 + . . . +X a =p. For 
instance, R 1234 =l!l!l!l!=l with a=4, and R 112223 =2!3!l!= 
2!3! with a=3. Some explanations of the general p-th order 
term in the last expansion in (8) are in order. It includes p+1 
types of terms ordered from the beginning to the end of the 
expression within the braces, with each type of terms except 
the first and last ones enclosed in square brackets. The c-th 
(l^c§p+l) type of terms is obtained by multiplying a 
(c-l)-th order term in the expansion of (<|> d k ) T by a (p-c+1)- 
th order term in the expansion of M(|> d k . For the c-th type of 
term, where 2^cSp, the p indices, i, j, . . . , and t are 
distributed in the subscripts of the two vectors pre- and 
post-multiplying M, whose numbers of indices in the sub- 
scripts are c-1 and p-c+1 respectively. For the first and last 
((p+l)-th) types of terms, all the p indices lie in the sub- 
scripts of the vectors post- and pre-multiplying M, respec- 
tively. The number of terms within each set of square 
brackets equals the number of different combinations of 
indices in the subscripts of the vectors pre- and post- 
multiplying M. When all the p indices, i,j, . . . , and t, have 
distinct values, due to symmetry of these vectors in their 
indices, terms of the c-th (l^c^p+1) type, involving dif- 
ferent permutations of the same indices in the subscripts of 
the vectors, are equal and combined, resulting in the scaling 
factor (c-l)!(p-c+l)! in front of the square brackets. Con- 
sequently, the indices in the second through p-th summations 
range from the values of their preceding summation indices 
to m. When any of the p indices, i, j, . . . , and t, have equal 
values, the corresponding terms in each type are given by 
those in the previous case divided by 



[0083] For instance, the 4 th order term of the form 
SG'oG^ in the expansion of (<|) d k ) T M(|) d k can be obtained 
from the expression for the p-th order term in (8): 



[0084] where p=4 and the four indices involved, i, j, 1, and 
o, are i=l and j=l=o=2. 

[0085] Equating the coefficients of the 8Gj (i=l, 2, . . . , m) 
terms in (8) and using symmetry of the mass matrix yields 

(<|> k ) T Mz (1)i k =0 (10) 

[0086] Equating the coefficients of the oGj&Gj terms and 
using symmetry of M and z (2)ij k yields 



[0087] for all i, j=l,2, . 
procedure, one obtains 



, m . Following a similar 



[0088] for i,j,l=l, 2, . . . , m. Equating the coefficients of 
the bG fiGj . . . bG e bG t terms of p-th order in (8) yield 



■■■ + (4/^ P -«-J + 

2 !(p - 2) ! [(4/m Z|p _ 2)1 ... „ + (4,,) 7 'm4_ 2)j . ... „ + 

■■■ + (4/^4-2)^..,] + ■■■ + 

(p-2)-2-[(4_ 2)j ... J ,)'M4^ + (4_ 2)J ,.. J ,) , M4, + 

-+(4-2w-/ M 4»] + 

- + ftU»-/"^]} 



4 ! ! # A/4 1222 + 3 ![(4/m4 222 + (4 2 ) 7 M4 1: 

(4l)2) 7 «4l22 + (42) 7M 4l22] + 

2 !2![(4 12 ) 7 M4 22 + (4 12 ) T M4 22 + 



[0089] for i,j, . . . ,t=l,2, . . . , m. The p-1 types of terms, 
enclosed in the p-1 sets of square brackets on the right-hand 
side of (13), are ordered from the beginning to the end of the 
e xpression within the braces, and their structures are readily 
observed. When p is odd, by symmetry of 
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M and z\ p)ij . . . , , the y - th{l • 



[0090] type of terms is identical to the (p-y)-th type, and 
the two types of terms can be combined. Similarly, when p 
is even, the 



[0097] is the coefficient of the u-th term and the number in 
the parentheses in its subscript corresponds to that of 



[0091] type of terms equals and can be combined with the 
(p-y)-th type of terms. In this case, however, there is a 
separate type, the 



[0098] Pre-multiplying (15) by (<j) k ) T and using (1), (2), 

and (16) yields 



[0092] of terms in the middle of the expression. 
[0093] Substituting (4)-(6) into (3) yields 



[0099] Substituting (16) into (10) and using (2) yields 



|a» + g A^.tfGi + g g ^ijSCiSCj + 
gg|^ C ,. aCj , Cl + ..j 

m|/ + g z^sg; + g g ^aG.aG,- + 



[0100] Pre-multiplying (15) by (<f) T , where 1 S V £N and 
v*k, and using (1), (2), and (16) yields 



[0101] By (16), (18) and (19) we have determined 



[0094] Equating the coefficients of the 6G ; (i=l,2, . . . , m) 
terms in (14) yields 



K : | + ' = I U-,, + I, i W 



[0095] Expanding z, ni k using normalized eigenvectors of 
(1) as basis vectors, one has 



[0102] Equating the coefficients of the terms in 

(14) yields 



t Oft L OA r 

-'l U-,;, + 1,1, + 1,1, +-'",:, »- A 



[0103] Expanding 
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[0104] using normalized eigenvectors of (1) as basis vec- 
tors, one has -continued 

2![Af : „,A^ I , + Mi 2ui M4 UJ +Af 2w M4 )1 ]3!A{ 3wl M/ 

N (21) 
Am = 2j *W [0112] Expanding 



[0105] where 



4* 



[0113] using normalized eigenvectors of (1) as basis v 
tors, one has 



[0106] is the coefficient of the u-th term and the number in ^ _ y ^ ^ ^ 

the parentheses in its subscript corresponds to that of 3,u ' ~i ' 



[0107] Pre-multiplying (20) by (<(> k ) T and using (1), (2), 
(10), and (21) yields 



[0114] where 



[0108] Substituting (21) into (11) and using (2) yields 



[0115] is the coefficient of the u-th term and the number ir 
the parentheses in its subscript corresponds to that of 



[0116] Pre-multiplying (25) by (<)> k ) T and using (1), (2), 
(10), and (26) yields 



[0109] Pre-multiplying (20) by (f ) T , where l^v^N and dK 

v*k, and using ( 1 ). (2), and (2 1 ) yields g^fa -^iv M 4m -^iMiw 



:ux -a- 

(\8K , dK 



(24) [0117] Substituting (26) into (12) and using (2) yields 



[0118] Pre-multiplying (25) by (f ) T , where 1 S v S N and 
[0110] By (21), (23) and (24) we have determined v *k, and using (1), (2), and (26) yields 



: vnr 



[0111] Equating the coefficients of the oG^-SOt terms 
(14) yields 



T dK dK dK 



3!A*M4 fel + 2![A a „-M4 2 , jl +A U)J Mzf 2)1 ,- + A^Mz^,,-,] + 



AfiyM^,,,. - AJdjM^j - Af 2)J ,.M z f„, - - A{ ;) , ; ,A^ 1)( ] 
(25) [0119] By (26), (28) and (29) we have determined 
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[0120] We proceed now to derive the perturbation solu- 
tions for the general p-th order terms in (5) and (6). Equating 
the coefficients of the bGfiG; . . . 8G s 8G t terms of order p in 
(14) yields 



[0126] The p-th order sensitivities of eigenvalues are 
obtained from (7) and (32). They depend on the eigenvalue 
and eigenvector sensitivities of orders up to p-2 and p-1 
respectively. Substituting (31) into (13) and using (2) \ ickls 



= p \X k M^ p)ij ... „ + (p - l)![Af lli M4_ 11J ,.. a + 
Afi, J M4_ 11 ,.. s , + -+Af 11 ,M4_ 11| ,,..J + 
2 !(p - 2) ![A* 21 , 7 m4_ 21 ,.. a + Af 2); ,M4_ 2)J ,.. a + 



[0121] Expanding 



[0122] using normalized eigenvectors of (1) as basis vi 
tors, one has 



(4>i/ M 4-2>;...» + ■■■ + {$2*7 M$ p _ 2)ll .. q \.. + 
(p - 2) !2![(4_ 2 „. ) T M4 2)ij + (4_ 2 , J „,,) 7 'M4 )il + 

- + (tr-w.jM^] + (p - \y\(Ap-m.J M A» + 

(Ap-m.J M A»i + ■■■ + (4-iw... s ) r M4„]} 



[0127] Pre -multiplying (30) by (<)> V ) T , where 1 Sv^N and 
v*k, and using (1), (2), and (31) yields 



" =1 2!(p-2)!(Af 2)i/ Aftf 2y ...„+, 



[0123] where 



^2)«M4_ 2w ... r ) - ... - (p - l)!(Af„_ 1)J ,.. s ,M4 )1 . + 

■4 i ),....>/ A '4i ), + ■■■+ -4-1 )<>■ ...j m 4i » 1 



[0124] is the coefficient of the u-th term and the number ir 
the parenthesis in its subscript corresponds to that of 



[0128] By (31), (33) and (34) we have determined 



[0129] The p-th order sensitivities of eigenvectors can be 
subsequently obtained from (7). They depend on the eigen- 
value and eigenvector sensitivities of orders up to p-1. 

[0125] Pre-multiplying (30) by (<^f and using (1), (10), 

(31) and orthonormality relations of eigenvectors yields t 0130 ] Equations (5) and (6) serve both the forward and 

inverse problems. In the former one determines the changes 
in the eigenparameters with changes in the stiffness param- 
k i T , QK (32) eters. Damage detection is treated as an inverse problem, in 

Apfi- si = y,^ 1 t(P - ^\gcAp-i)j- " + which one identifies iteratively the changes in the stiffness 

gK gf , parameters from a selected set of measured eigenparameters 

^-4-ni n + ■ ■ ■ + -g^Ap-m - «j~ of the damaged structure: k d k (k=l, 2, . . . , %) and <|) d k (k=l, 

1 2, . . . , n^,), where l = n x , n^N. Here X d k and (|) d k are 

(p-l)!(Af 1)i M4-Dj w 4-ni « + - + assumed to be the perfect eigenparameters; simulated noise 

, 4 ,, k , ... „..,,* ,,_t is included in the measured eigenparameters in the beam and 

frame examples that follow. Often we choose a set of n 

i:,Ur_ : Jf + ...+ I,, ,,.!/- _! ,i measured eigenparameter pairs to detect damage, i.e., 

. _2)nux k M?* + %=n+=n. Let the number of the measured degrees of free- 

P ^ 2) '- s: ™ dom of (|) d k be N m ; N m =N and N m <N when we have 

-4-2)j---« M 4>i! + --- +A (p-2>!/-- r Mz a)«)] complete and incomplete eigenvector measurements, 

respectively. With reduced measurements the unmeasured 
degrees of freedom of <|) d k is estimated first using a modified 
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eigenvector expansion method (see the beam and frame 
examples below) and <|) d k is mass-normalized subsequently. 
Only the component equations corresponding to the mea- 
sured degrees of freedom of <)) d k are used in (6). The system 
equations in Eqs. (5) and (6) involves n x +n^N m scalar 
equations with m unknowns, which are in general determi- 
nate if n x +n (|) N m =m, under-determined if n^+n (( ,N m <m, and 
over-determined if %+n^N^m. In the first iteration, X k and 
<|) k in (5) and (6) correspond to the eigenparameters of the 
structure with the initial stiffness parameters G^. With the 
perturbation terms determined as shown above, the changes 
in the stiffness parameters 8G ; (i:i , where the number in the 
superscript denotes the iteration number, can be solved from 
(5) and (6) using an optimization method discussed below. 
The estimated stiffness parameters of the damaged structure 
are updated by G ; (1) =G ; (0) +8G ; (1) . With the updated stiffness 
parameters, the eigenparameters, X k and <() k , in (5) and (6) are 
re-calculated from the eigenvalue problem (1) and the 
perturbation terms are re -evaluated. One subsequently finds 
6G ; (2) using the same method as that in the first il 
This process is continued until the terminal 
|5Gj 1 |<e, where L is the last iteration number and e is some 
small constant, is satisfied for all i=l,2, . . . , m. Note that 
after the w-th (l^w<L) iteration, we set G ; (w) to G hi if 
G ; (w:> >G h; , and to zero or some small stiffness value e G if 
Gj (w) <0. A flowchart for the iterative algorithm is shown in 
FIG. IB. When a single iteration is used, the iterative 
method becomes a direct method. 

[0131] FIG. IB shows steps included in a method for 
detecting structural damage in accordance with one embodi- 
ment of the present invention. The method includes as an 
initial step measuring one or more eigenparameters, k d k , <|> d k 
(Block I). These eigenparameters are then compared with 
estimated eigenparameters associated with the stiffness 
parameters, G^ (Block 2). The differences between the 
measured and estimated eigenparameters are then used by 
the perturbation method to establish system equations (5) 
and (d) (Block 3). The perturbation method may be a first or 
higher order multiple-parameter perturbation procedure. 

[0132] Next, an optimization method may be used to find 
the changes in the stiffness parameters oG ; lw) (Block 4). 
These values are then compared to the predetermined value 
e (Block 5). and if the absolute values are less than e the 
stiffness parameters identified are set as G i fw " 1) (Block 6). 
[0133] If the absolute values are not all less than e, the 
process proceeds along an iterative path where the stiffness 
parameters are first updated (Block 7). The stiffness param- 
eters are then bounded between 0 or e G and G hi (Block 8), 
and eigenparameters associated with the updated stiffness 
parameters are calculated (Block 9). The comparison of 
Block 2 is then performed based on these calculated, or 
estimated, eigenparameters. 



[0134] Opt 



Methods 



[0135] Neglecting the residuals of order p+1 in (5) and (6) 
yields a system of polynomial equations of order p. When 
n x+ n $ N m = m - 8G i (w) (i=l,2, . . . , m) at the w-th iteration can 
be solved from the resulting equations. There are generally 
an infinite number of solutions when n x +n l|) N m <m, a unique 
solution when n x +n (|) N m =m and p=l, and a finite number of 
solutions when n x +n (|) N m =m and p>l. When n^+n^N^m, 
one generally cannot find 6G ; (w) to satisfy all the equations, 
and an optimization method can be used to find 8Gj <w) which 



objective function related to the 6 
satisfying these equations. We use here the s; 
e x k and e^, to denote the errors in satisfying the system 
equations in (5) and (6) respectively. Consider the objective 
function 



[0136] where W x k (k=l,2, . . . , n,) and W^ k (k=l, 2, . . . 
, n^) are the weighting factors, and J is a function of 8Gj (w) 
when one substitutes the expressions for e x k and e^ in (5) 
and (6) into (35). When the first-order perturbations are 
retained in (5) and (6), the least-squares method can be used 
to determine SG^ which minimize (35) with W x k =W„, k =l. 
Other weighting factors can be handled bv dividing (5) and 
(6) by 



-== ana -== 



[0137] respectively. The method determines essentially 
the generalized inverse of the resulting system matrix, and 
is also known as the generalized inverse method. When the 
perturbations up to the p-th (p=l) order are included in (5) 
and (6), the gradient and quasi-Newton methods [25] can be 
used to determine 6G ; lw) iteratively. Unlike the generalized 
inverse method, the methods are applicable when other 
objective functions are defined. While the optimization 
methods are introduced for over-determined systems, they 
can be used when n^+n^N^m, in which case J=0 (i.e., 
e^e^O) when the optimal solutions are reached. 

[0138] Gradient Method 

[0139] To minimize the objective function in (35) at the 
w-th iteration, one can use the algorithm 



[0141] a b >0 is the step size, and f b _ 1 equals the gradient 



US 2005/0072234 Al 



11 



Apr. 7, 2005 



[0142] associated with 



[0143] Note that the subscript b (b & 1) in all the variables 
in (36) denotes the number of nested iterations. The initial 
values used are 



the quasi-Newton methods, including the step size search 
procedure as outlined below, is shown in FIG. 2. 

[0150] Step Size Search Procedure 

[0151] The optimal step size for the gradient and quasi- 
Newton methods is determined in each nested iteration to 
minimize the function 



A* 50 !*-!) -<*a/a-i ) = F ( a b) 



[0144] The nested iteration is terminated when a b ||g b _ 
i|U<Y> where H'l^ is the infinity norm and y is some small 
constant, or the number of nested iterations exceeds an 
acceptable number, D. 

[0145] Quasi-Newton Methods 

[0146] Due to its successive linear approximations to the 
objective function, the gradient method can progress slowly 
when approaching a stationary point. The quasi-Newton 
methods can provide a remedy to the problem by using 
essentially quadratic approximations to the objective func- 
tion near the stationary point. The iteration scheme of these 
methods is given by (36) with f b .i=B b . 1 g b . 1 , where B^ is an 
approximation to the inverse of the Hessian matrix used at 
the b-th nested iteration, and the other variables the same as 
those previously discussed. Initially, we set 



<SGS$ = 0 



[0147] and B 0 =I, the identity matrix. The matrix B b is 
updated using either the DFP formula 

-,)f,;, , nor,;, , / < 37 ) 

lB b -i (gb - gb-l )] Wb-i (gb - gb-l)f 
(gb-gb-lfBb-iigb-gb-O 



[0152] with respect to a b . The search procedure is divided 
into two phases: an initial search to bracket the optimum a* b 
and a golden section search to locate a* b within the bracket, 
as shown in FIG. 2. 

[0153] Initial Bracketing. Choose the starting point Xj=0 
and an initial increment A>0. Let x 2 =x 1 +&, F 1 =F(x 1 ), and 
F 2 =F(x 2 ). Since for the gradient and quasi-Newton methods, 
f 0 =g 0 and it is along a descent direction of J when A is 
sufficiently small, one has F 2 <F 1 . Rename 2A as A, and let 
x 3 =x 2 +A and F 3 =F(x 3 ). If F 3 >F 2 , stop and a* b is contained 
in the interval (x 1; x 3 ). Otherwise, rename x 2 as x 1 and x 3 as 
x 2 , then F 2 becomes F 1 and F 3 becomes F 2 . Rename 2A as A, 
and let x 3 =x 2 +A and F 3 =F(x 3 ). Compare F 3 and F 2 , and 
repeat the above procedure if F 3 <F 2 until F 3 >F 2 , with the 
final interval (x 1; x 3 ) containing a* b . 

[0154] Golden Section Search. If |x 3 -x 2 |>|x 2 -x 1 |, define a 
new point: 

x 4 =x 2+ 0.382(x 3 -x 2 ) (39) 

[0155] Otherwise, rename x 1 as x 3 and x 3 as x 1; and then 
define x 4 using (40). Let F 4 =F(x 4 ). If F 2 <F 4 , rename x 4 as x 3 , 
then I , becomes l\. Otherwise, rename x, as x, and x 4 as x 2 , 
then F 2 becomes F 1 and F 4 becomes F 2 . Compare |x 3 -x 2 | and 
|x 2 -Xj|, and repeat the above procedure until |x 3 -xj|f b _ 
i|U <eol , where e tt is some small constant. Then choose 




[0148] or the BFGS formula 

B b = Bb-i + \i + = — '—^ ~^]\ 

(scft'-sd^/igb-gb-i) 

(6G t $-6G™ l) )(gb-gb- l ) T Bb- l + 
B b -l(gb -gb-l)(gb -gb-lf 
((S^-SG^vfigb-gb-O) 

[0149] The nested iteration is terminated when aJIB^g,,. 
i , <y or the number of iterations exceeds D. A flowchart for 



MASS-SPRING SYSTEM EXAMPLE 

[0156] The algorithm discussed above is used to identify 
the stiffness parameters of a N-degree-of-freedom system 
consisting of a serial chain of masses and springs, such as the 
system shown in FIG. 3. Let the masses of the system be 
M;=l kg (i=l, 2, . . . , N), and the stiffnesses of the 
undamaged springs be G h; =lN/m (i=l, 2, . . . , m), where 
m=N+l. The system is said to have a small, medium and 
large level of damage if the maximum reduction in the 
stiffnesses is within 30%, between 30 and 70%, and over 
70%, respectively. The mass matrix M is an NxN identity 
matrix, and the stiffness matrix K is a banded matrix with 
entries K^G^G^ (i=l,2, . . . , N), K^-K^—G^ 
(i=l, 2, . . . , N-l), and all other entries equal to zero. The 
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[0157] have a u 
respectively, and zt 
of the matrices 



Since vanishing stiffness in any spring other than the two 
end ones in FIG. 1 can result in two decoupled subsystems, 
when G ; (w) <0 we set G ; (w) to e G =0.1N/m in the first two 
iterations and to O.lN/m in the remaining iterations for all 
the cases considered here. A relatively large value is 
assigned to e Q in the initial iterations to avoid close eigen- 
values in the mass-spring system and small denominators in 
(19). This improves convergence especially when a small 
number of eigenparameter pairs are used. A smaller value is 
used for e G in the later iterations to improve the accuracy of 
stiffness estimation when there is a large level of stiffness 
reduction. Using the first-order perturbations and different 
numbers of eigenparameter pairs, the maximum errors in 
estimating the stiffnesses of the damaged system at the w-th 
iteration, defined by 



[0158] We look at a forward problem first with N=3 and 
m=4. The stiffnesses of the damaged system are G dl =G d3 = 
lN/m, G d2 =0.3N/m and G d4 =0N/m. The undamaged system 
is considered as the unperturbed system and the damaged 
system as the perturbed system. Based on the eigenparam- 
eters of the undamaged system, the eigensolutions of the 
damaged system are obtained using the first-, second-, and 
third-order perturbations, as shown in Table 1 below. The 
results show that even with the large changes in stiffness, the 
third-order perturbation solutions compare favorably with 
the exact solutions for the damaged system. The higher- 
order perturbation solutions can be used for large order 
systems when their direct eigensolutions become costly. 



[0160] are shown in FIGS. 4A and 4B for all the itera- 
tions. In FIG. 4A p=l and n- 1,2,3, and in FIG. 4B p=l and 
n=4,5, ... ,9. When n=l, the error decreases slow ly, though 
monotonically, and there is an estimation error of 1.5% at the 
end of iteration. While the errors can increase with the 
iteration number before approaching zero for n=3, they 
decrease monotonically for n=2 and n^4. All the stiffnesses 
are exactly identified at the end of iteration when n^2. Note 
that the number of the system equations equals and exceeds 
the number of unknowns when n=l and n^2, respectively. 
Since the system equations are linear, they have a unique 



| n.lf,Sih-| 
J ii. 05733 
1 11.73528 I 



TABLE 1 



n.05450 \ 
0.10607 1 
-0.45962 J 



o. 34478 -i 
-0.74053 I 
0.61507 J 



j 0.5 | 

J -0.70711 \ 
{ 0.50000 J 



{ 0.296351 
-0.74132 1 
0.62482 J 



{ 0.26445"! 
-0.74875 I 
0.62362 J 



[0159] Consider now the damage detection problem with 
N=9, m=10, G d5 =0.5N/m, G d8 =0.7N/m, G dlo =0.8N/m and 
G d; =lN/m (i=l,2,3,4,6,7,9). We set W^W^-l, e=0.001, 
y=10~ 10 , a ;-1 for all i, and D=500; the actual numbers of 
nested iterations in all the cases are much smaller than D. 



solution when n=l, and J has a unique minimum when n^2. 
With the small y the gradient method and the quasi-Newton 
methods using the DFP and BFGS formulas yield exactly the 
same results as the generalized inverse method (not shown 
here). Because the generalized inverse method does not 
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involve any nested iteration, it is the most efficient one 
among the four methods. While not shown here, the results 
indicated that that the quasi-Newton methods converge 
faster than the gradient method and the BFGS method has 
the similar performance to the DFP method. In what follows 
the BFGS method will be used with the higher-order per- 
turbations. With the second-order perturbations the errors 
shown in FIG. 4C decrease monotonically for all n (n=l,2, 
. . . , 9). The errors at w=l in the expanded view decrease in 
the order n=3,2,4,9,7,8, with the lines for n=5 and 7 virtually 
indistinguishable. In FIG. 4(A-C) the following symbols are 
used for different n: -—407 n=l; n=2; * n=3; 

— n=4; — n=5; — n=6; — n=7 * n=8; - n=9. While 
the use of the second-order perturbations improves the 
accuracy of stiffness estimation in each iteration and reduces 
the number of iterations, it takes a much longer time to 
compute the higher-order perturbations and the associated 
optimal solutions. 

[0161] When only the first few eigenvalues are used, for 
instance, n x =5 and n^=0, the stiffnesses identified with the 
first-order perturbations, 

0 926 ( 0 976 °0 9 875)°' N/m °' 8M ' 0 - 6 "' 0 - 6 "' ° MA ' 41 

[0162] where the number in the superscript denotes the 
last iteration number, correspond to those of a different 
system with the same eigenvalues for the first five modes as 
the damaged system. The same stiffnesses are identified with 
the second-order perturbations. Similarly, when the first 
eigenvector is used, i.e., ^=1 and n x =0, the stiffnesses 
identified with the first-order perturbations are those of a 
different system with the same eigenvector for the first mode 
as the damaged system: 

G< 9 >=(0.989, 0.991. 0.995, 1. 0.520. 0.829, 0.940, 

0.668, 0.959, 0.768) T N/m (42) 

[0163] With the second-order perturbations the stiffnesses 
of the damaged system are identified. The stiffnesses iden- 
tified are not unique because the system equations in each 
iteration are under-determined. The solution given by the 
generalized inverse method here in each iteration is the 
minimum norm solution. Increasing the number of eigen- 
parameters used can avoid this problem. 

[0164] If the system has a large level of damage, i.e., 
G d5 =0.3 N/m, G dlo =0.1 N/m, with the other parameters 
unchanged, the stiffnesses of the damaged system are iden- 
tified with the first-order perturbations after 55 iterations 
when n=l and 6 iterations when n=2 and 3 as shown in FIG. 
5A. For n=4,5,6,7, the errors decrease monotonically and the 
number of iterations is reduced slightly, as shown in FIG. 
5B, which has an expanded view near w=l. With the 
second-order perturbations the errors shown in FIG. 5C 
decrease monotonically for n=l,2, . . . , 5, and the number 
of iteration for n=l is reduced from 55 in FIG. 5A to 4. The 
errors at w=l in the expanded view in FIG. 5C decrease in 
the order n=3,2,4,5, with the lines for n=2 and 4 virtually 
indistinguishable. The symbols used in FIG. 5A-C for n=l, 
2, . . . , 7 are the same as those in FIG. 4A-C. 

[0165] Finally, consider a large order system with a large 
level of damage: N=39, m=40, G dl2 =0.7N/m, G dl9 =G d37 = 
0.1 N/m, G d28 =0.8 N/m, G dl =l N/m (1^40 and i* 12,19, 
28,37), and the other parameters are the same as those in the 
previous example. With the first-order perturbations the 
exact stiffnesses arc identified after 57 iterations when n=l, 



as shown FIG. 6A. Using a larger number of eigenparameter 
pairs (n=2,3,4,5) significantly reduces the number of itera- 
tions, as shown in FIG. 6B. The errors at w=2 in the 
expanded view in FIG. 6B decrease in the order n=4,2,3,5. 
With the second-order perturbations the exact stiffnesses are 
identified after 6 iterations when n=l and 3 iterations when 
n=2, as shown in FIG. 6C, which has an expanded view for 
1£ W S4. The symbols used in FIG. 6(A-C) for n=l, 2, . . . 
, 5 are the same as those in FIGS. 4A-C and 5A-C. 

FIXED-FIXED BEAM EXAMPLE 

[0166] The algorithm discussed above may be applied to 
detecting structural damage in an aluminum beam with fixed 
boundaries. The beam of length L t =0.7 m, width W=0.0254 
m, and thickness H=0.0031 m has an area moment of inertia 



/ = — WH 5 = 6.3058 x 10 -11 m 4 



[0167] and a mass density p=2715 kg/m 3 . The finite 
element method is used to model its transverse vibration. 
The beam is divided into N e elements, as shown in FIG. 7, 
with the length of each element being 




[0168] There are N e +1 nodes. With V ; and 6; denoting the 
translational and rotational displacements at node i (i=l, 2, 
. . . , N e +1), the displacement vector of the i-th (i=l,2, . . . 
, N e ) element is [Vi,e;,V ;+1 ,6 i+1 ] T . The Young's modulus is 
assumed to be constant over each beam element and that of 
the i-th element is denoted by G ; . The Young's modulus of 
the undamaged beam is G h =69xl0 9 N/m 2 . Hence G h; =G h for 
i=l,2, . . . , m, where m=N e . Small to large levels of damage 
correspond to reductions of the moduli defined for the 
Mass-Spring Example. The mass and stiffness matrices of 
the i-th beam element are 





156 


-224 54 


134 


PWHI, 


-224 


4% -134 


-■ill 




54 


-134 156 


22/, 




134 


-31] 224 






[0169] Using the standard assembly process yields the 
2(N e +l)x2(N e +l) global mass and stiffness matrices. Con- 
straining the translational and rotational displacements of 
the two nodes at the boundaries to zero yields the NxN M 
and K matrices, where N=2(N e -l) is the degrees of freedom 
of the system. The displacement vector of the system, 
involving the displacements of the 2 nd through N e -th node, 
is [V 2 ,e 2 ,V 3 ,e 3 , . . . ,V Ne , e N J T . The matrix 
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dK 



[0170] (i=l, 2, . . . , m) can be obtained from K by setting 
G ; =l and G 1= . . . =G ; _ 1 =G ;+1 = . . . =G N =0. The parameters 
W x k , W^, e, y, a ; , and D are set to the same values as those 
previously discussed, and e G is set to 0.15 G h in the first two 
iterations and to 0.05 G h in the remaining iterations. The 
first-order perturbations are used unless indicated otherwise. 

[0171] Consider first the cases with N e =m=10 and N-18. 
When the system has a medium level of damage: 

G d =(l, 1, 1, 1, 0.5, 1,1, 0.7, 1, 0.8) T xG„ (44) 

[0172] the stiffness parameters of the damaged system are 
identified after 6 iterations with n=l. When the system has 
a large level of damage: 

G d =(l, 1, 1, 1, 0.3, 1, 1, 0.7, 1, 0.1) T xG h (45) 
[0173] the stiffness parameters of the damaged system are 
identified after 7 iterations with n=l. Consider next the cases 
with N e =20, 40 and 80. For the systems with medium and 
large levels of damage, the stiffness parameters of the first 10 
elements are given by (44) and (45), respectively, and those 
of the remaining elements are G h . In all the cases the 
stiffness parameters of the damaged systems are identified 
within 10 iterations when n=l. The numbers of iterations are 
reduced slightly when the second-order perturbations are 
used. Note that the system equations are over-determined 
when n=l. 

[0174] When only the translational degrees of freedom of 
an eigenvector are measured, a modified eigenvector expan- 
sion method is used to estimate the unmeasured rotational 
degrees of freedom. To this end, (() d k is partitioned in the 
form <|> d k =[(<|> dm k ) T , (^d^Y* where (|> drn k and ()> du k are the 
measured and unmeasured degrees of freedom of <|> d k , 
respectively. Similarly, <|> k in (6) is partitioned in the form 
<l )k =[( < l ) m k ) T ' ( ( l ) u k ) T ] T ' where (|) m k and (|> u k correspond to the 
measured and unmeasured components of (() d k , respectively. 
Since <|) dm k , <|) m k and (|) u k are known in each iteration, <|> du k is 
estimated from <|> du k =k[(i|> m k ) + (|) dm k ](|) u k , where the super- 
script+denotes generalized inverse. Once the rotational 
degrees of freedom of (|) d k are determined, (() d k and (j) k are 
converted to their original forms and (() d k is mass-normal- 
ized. Only the component equations corresponding to the 
measured degrees of freedom of (|) d k are used in (6) and the 
system equations are determinate when n=l. The exact 
stiffness parameters of the damaged systems considered 
above can be identified. For the 10- and 20-element beams 
with the medium levels of damage, the stiffness parameters 
are identified after 6 and 24 iterations, respectively, when 
n x =l and ^=2. For the 10- and 20-element beams with the 
large levels of damage, the stiffness parameters are identified 
after 9 iterations when n x =l and ^=3 and 10 iterations when 
n=2, respectively. 

[0175] Finally, the effects of measurement noise on the 
performance of the algorithm are evaluated for the 10-ele- 
ment beam with the large level of damage. Simulated noise 
is included in the measured eigenparameters: 

X d k =X* d k +v^ k X* d k , W^rZ+v-^Vd 11 (46) 



[0176] where X* d k and (|)* d k are the k-th perfect eigenvalue 
and eigenvector, respectively, R x k is a uniformly distributed 
random variable in the interval [-1,1], R^ k is a diagonal 
matrix whose diagonal entries are independently, uniformly 
distributed random variables in the interval [-1,1], and 
ve[0,l] is the noise level. Note that R x k and are gener- 
ated for each measured mode. Each random parameter is 
generated 10 times and the average is used. Three different 
noise levels are considered: v=5'i, 10' < and 20' , . When all 
the degrees of freedom of an eigenvector are measured, the 
stiffness parameters identified with n=l, 2, and 3 are shown 
in FIGS. 8A, 8B, and 8C, respectively. The following 
symbols are used for different noise levels: v=5%; Q , 
v=10%; B, v=20%. When only the translational degrees of 
freedom of an eigenvector are measured, the eigenvector 
expansion method described above is used and the stiffness 
parameters identified with n=2 and n=3 are shown in FIGS. 
9A and 9B, respectively. The same symbols as those in FIG. 
8A-C are used in FIG. 9A-B for different noise levels. The 
stiffness parameters corresponding to v=0 in FIGS. 8 and 9 
are the exact values. It is seen that in the presence of noise, 
the stiffness parameters can be accurately identified with an 
increased number of measured eigenparameters. 

[0177] Space Frame Example 

[0178] The damage detection method can be applied to 
more complex structures, such as. for example, the modular, 
four bay space frame shown in FIG. 10. In this example and 
all the following examples the first order perturbation 
approach along with the generalized inverse method is used. 
The four-bay space frame example shown in FIG. 10 is 
made of extruded aluminum with 48 members, and is 8' tall, 
1.46' wide and 1.8' deep. The horizontal and diagonal 
members have the same cross-sectional dimensions (25.4 
mmxl2.7 mm) and the vertical members are "L"-angles 
with cross-sectional dimensions 50.8 mmx50.8 mmx6.35 
mm (thickness). All the members are connected through 
bolted joints. The frame is assumed to be fixed to the ground. 
The finite element method is used to model the 3-dimen- 
sional vibration of the frame, with each member modeled 
with four 12-degrees-of-freedom beam elements. The total 
degrees of freedom are 960 when the boundary conditions 
are applied. The Young's modulus is assumed to be constant 
over each member and that of the i-th (1 ^i£m=48) member 
is denoted by G ; . The Young's modulus of the undamaged 
member is G h =69xl0 9 N/m 2 and G h; =G h . A vertical (mem- 
ber 1) and a diagonal (member 48) member in the first bay 
are assumed to have 70% and 30% reductions in Young's 
modulus, respectively, and all the other members are 
assumed to be undamaged. The first two vibration modes, 
corresponding to the bending of the frame along the x and 
y directions, respectively, are used to detect damage. The 
translational degrees of freedom of the 16 nodes, denoted by 
"S" in FIG. 10, along the x and y directions are assumed to 
be measured with 10% measurement noise. All the other 
degrees of freedom are estimated using the eigenvector 
expansion method discussed above and the measured modes 
are mass-normalized. The Young's moduli of all the mem- 
bers of the damaged truss are identified as shown in FIG. 10, 
with the maximum estimation error less than 7%. Note that 
the truss has closely spaced vibration modes. With the 
modes arranged in the order of increasing frequencies in 
each iteration, the method has effectively handled mode 
switching that has occurred in the damage detection process. 
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[0179] In summary, the damage detection method identi- 
fies stiffness parameters in structures, which have a small, 
medium, and large level of damage if the maximum reduc- 
tion in the stiffnesses is within 30%, between 30 and 70%, 
and over 70%, respectively: A large level of damage is 
studied in many examples because this poses the most 
challenging case, with sever mismatch between the eigen- 
parameters of the damaged and undamaged structures. 

[0180] The damage detection method as embodied and 
broadly described herein can be applied to structures that 
can be modeled with beam elements. A beam element is an 
element that has one dimension that is much longer than the 
other two. This element is very good at modeling "I"-beams, 
rectangular beams, circular beams, "L"-angles, "C"-chan- 
nels, pipes, and beams with varying cross sections. Struc- 
tures that can be modeled with this element include, but are 
not limited to, lightning masts, light poles, traffic control 
poles, pillar type supports, bridges, pipelines, steel building 
frameworks, television, radio, and cellular towers, space 
structures, cranes, pipelines, railway tracks, and vehicle 
frames. Structures that can be modeled as beam elements are 
used simply for case of discussion, and it is well understood 
that the damage detection method discussed above may be 
applied to structures that can be modeled by other elements 
using the finite element method or modeled by using other 
methods. 

[0181] Damage Detection Using Changes of Natural Fre- 
quencies: Simulation and Experimental Validation 

[0 1 82] I 'or structures such as beams and lightning masts in 
electric substations, using only the changes in the natural 
frequencies can relatively accurately detect the location(s) 
and extent of damage, even though the system equations are 
severely underdetermined in each iteration. This is an inter- 
esting finding as it is much easier to measure the natural 
frequencies than the mode shapes, and demonstrates the 
effectiveness of the iterative algorithm. Extensive numerical 
simulations on beams and lightning masts confirmed this 
finding. Experiments on the beam test specimens with 
different damage scenarios and a lightning mast in an 
electric substation validated the simulation results. The 
beam test specimens and the lightning mast are used as 
examples for demonstration purposes, and the method can 
be applied to other structures. Note that unlike the beam 
example shown earlier, where the Euler-Bernoulli beam 
finite element model is used, the Timoshenko beam finite 
element model is used in all the examples here. The Timosh- 
enko beam theory is found to be more accurate in predicting 
the natural frequencies of the lightning masts and circular 
beams than the Euler-Bernoulli beam theory. 

[0183] For a cantilever beam, simulation results show that 
the damage located at a position within 0-35% and 50-95' ! 
of the length of the beam from the cantilevered end can be 
easily detected with less than 5 measured natural frequen- 
cies, and the damage located at a position within 35-50% of 
the length of the beam from the cantilevered end and at a 
position within 5% of the length of the beam from the free 
end can be relatively accurately detected with 10-15 mea- 
sured natural frequencies. 



[0184] Numerical and Experimental Verification 

[0185] Cantilever Aluminum Beams 

[0186] Experimental damage detection results for four 
different scenarios are shown first, followed by various 
simulation results. 

[0187] Scenario 1: Evenly-Distributed Damage Machined 
from the Top and the Bottom Surfaces of the Beam Test 
Specimen. 

[01 88] The aluminum beam test specimen shown in FIG. 
12 is 45 cm long by 2.54 cm wide by 0.635 cm thick. It is 
divided into 40 elements (each element has a length of 1 . 1 25 
cm). The beam has a section (from approximately 10 cm to 
15 cm from the cantilevered end) of 5 cm long and 7.62E-4 
m thick machined both from the top and the bottom surfaces 
of the beam. This corresponds to 56% of damage (or 
reduction of bending stiffness EI) along the length of five 
elements (from the 9 th to the 13 rd element). Using the 
changes of the first 2 to 5 measured natural frequencies, 
damage is detected within 7 elements using 2 or 5 measured 
frequencies (from the 7 th to the 13 rd element with 5 mea- 
sured frequencies and from the 5 th to the 11 th element with 
2 measured frequencies) and within 8 elements using 3 or 4 
measured natural frequencies (both from the 6 th to the 13 rd 
element), as shown in FIG. 13. With 5 measured frequencies 
the average damage of the 7 elements in FIG. 13, from 6.75 
cm to 14.625 cm, is 46%. Note that the elements with 
damage less than 10% are not accounted for here. The extent 
of damage detected is slightly lower than the actual extent 
because the predicted damage occurs at 2 mote elements (the 
7 th and 8 th elements) than the actual one. The error results 
from the solution of the severely underdetermined system 
equations (5 equations with 80 unknowns). 

[0189] Scenario 2: The Same Aluminum Beam Test Speci- 
men as in Scenario 1 Clamped at the Other End. 

[0190] The same beam was tested in the clamped-free 
configuration with the clamped end reversed, which placed 
the damage from 25 cm to 30 cm (from the 23 rd to the 27 th 
element) from the cantilevered end. The damage detection 
results with 3 to 5 measured frequencies are shown in FIG. 
14. With 5 measured frequencies the average damage 
between 23.625 cm to 32.625 cm (from the 22 nd to the 29 th 
element) is 40%. Again the elements with damage less than 
10% are not included. 

[0191] Scenario 3: Undamaged Cantilever Aluminum 
Beam Test Specimen with the Same Dimensions as Above. 

[0192] An undamaged aluminum beam test specimen was 
clamped at one end with the same configuration as shown in 
FIG. 12. With the first 3 to 4 measured frequencies the 
maximum error for the estimated bending stiffnesses of all 
the elements is within 10%, as shown in FIG. 15. 

[0193] Scenario 4: A Cut of Small Width on a Cantilever 
Aluminum Beam Test Specimen, Shown in FIG. 16 with the 
Same Dimensions as Above. 

[0194] The beam shown below has a cut that is 0.4191 cm 
deep and 0.1016 cm wide, which corresponds to a 96% 
reduction in bending stiffness at the cut. The beam is divided 
into 45 elements and the cut is located in the middle of the 
23 rd element. With 3 to 5 measured natural frequencies, the 
damage is detected at several elements surrounding the 23 rd 
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element, as shown in FIG. 17. It was found there was a 
roughly 50%, 40%, and 60% reduction in stiffness at the 
22'"', 23"', and 24 lh element, respectively. The estimated 
damage extent at the above elements is lower than that at the 
cut because the damage is distributed along these elements. 
[0195] To examine the effectiveness and robustness of the 
damage detection algorithm, various simulations with dif- 
ferent damage scenarios were carried out. In this way, we 
can gain more insight concerning the accuracy of the finite 
element model, convergence of the estimated bending stiff- 
nesses of all the elements of the beam with the increased 
numbers of measured natural frequencies and/or mode 
shapes, and region of the beam within which the damage can 
be detected with few measured frequencies. With the can- 
tilever aluminum beam divided into 40 elements, the fol- 
lowing two simulations have the similar damage location 
and extent to those in Scenarios 1 and 2 in the experiments: 
[0196] Simulation 1: Uniform Damage Between 9.0 cm 
and 15.75 cm from the Cantilevered End of the Beam with 
the Same Dimensions as Discussed Above. 

[0197] Simulation 2: Uniform Damage Between 29.25 cm 
and 36 cm from the Cantilevered End of the Beam with the 
Same Dimensions as Discussed Above. 
[0198] Simulation results in FIG. 18 and FIG. 19 show 
that the damage can be relatively accurately detected with 
the first 3 to 5 measured frequencies, and can be accurately 
detected with the first 10 measured frequencies. 
[0199] With the cantilever aluminum beam divided into 
the same number of elements, two more simulations are 
presented here: one for a multiple damage scenario — 70% of 
damage at the 3 rd element and 30% of damage at the 20 th 
element, and the other for a 50% of uniform damage from 
the 16 th to the 18 th element (i.e., the damage is located at a 
position within 35-50% of the length of the beam from the 
cantilevered end). 

[0200] Simulation 3: Multiple Damage of the Beam with 
the Same Dimensions as Discussed Above: 70% of Damage 
at the 3 rd Element and 30% of Damage at the 20 lh Element. 
[0201] Simulation results in FIG. 20 show that the mul- 
tiple damage can be relatively accurately detected with the 
first 3 to 5 measured frequencies, and accurately detected 
with the first 10 measured frequencies. 
[0202] Simulation 4: Uniform Damage from the 16 lh to the 
18 th Element of the Beam with the Same Dimensions as 
Discussed Above. 

[0203] Simulation results in FIG. 21 show that the dam- 
age within 35-50% of the length of the beam from the 
cantilever end can be accurately detected with 10-15 mea- 
sured frequencies. 

[0204] Lightning Masts 

[0205] The lighting mast shown in FIG. 22 has two 

sections of constant cross-sections and an eccentric spike. 
The lengths of the lower and upper sections are 6.89 m and 
6.83 m, respectively, and that of the spike above the upper 
section is 1.4375 m. Both mode shapes and natural frequen- 
cies were measured for this mast. The mode shapes were 
measured by using a laser Doppler vibrometer. The natural 
frequencies were measured using both the laser Doppler 
vibrometer and an accelerometer. It takes about 0.5-1.5 days 



the mode shapes and about 30 minutes to 
the natural frequencies. It would be even more 
difficult to measure the mode shapes for some of the taller 
lightning masts, which are 100 and 130 feet tall. A finite 
element model was made using OpenFEM. Once the model 
was completed the measured and calculated natural frequen- 
cies and mode shapes were compared. The mast was 
expected to be undamaged, since by inspection the measured 
and calculated natural frequencies matched (Table 2). It was 
found that the first mode shape measurement is affected by 
wind and high modes are less affected by the wind. The 
measured frequencies are not affected by the wind and can 
be used for damage detection. 

TABLE 2 

Comparison of mcasiucj and calculated (lioni Ihc (mile etcmcnl 



-2.4104 
-0.8418 
-6.1104 
0.31327 

29.26 30.36 -3.7385 

45.32 47.85 -5.5831 

47.75 50.87 -6.5191 

54.10 55.81 -3.1439 

67.45 68.83 -2.0476 

74.70 78.99 -5.7456 



[0206] Damage detection was then performed using only 
the first 4 to 7 measured natural frequencies, as shown in 
Table 2. The mast is modeled by 40 elements; the lower 
section has 18 elements, the upper section has 15 elements, 
and the spike has 7 elements. The experimental damage 
detection results are shown in FIG. 23. The mast is modeled 
by 40 elements with 20 elements in each of the two sections. 
It appears that the elements around the joints near the middle 
of the mast (8.2 m in FIG. 23) and at the free end of the mast 
(13.72 m in FIG. 23) had some damage. This is most likely 
due to the fact that the joints were modeled in the finite 
element model as being infinitely stiff, which is almost never 
the case. Also, the spike is actually connected to the free end 
of the mast at two points. In the finite element model here 
it is assumed that the spike is attached to the mast along the 
whole line of contact that has a length of 0.34 m, which 
makes the model a little stiffer. Similar damage detection 
results were obtained with different numbers of measured 
frequencies. 

[0207] Simulations for the same lightning mast as shown 
in FIG. 22 with different damage scenarios are carried out. 
Shown in FIG. 24 are the simulation results for the mast 
with 70% damage between 0.76 m and 1.15 m from the 
ground. Shown in FIG. 25 are the simulation results for the 
mast with 50'r damage between 6.89 m and 7.35 m from the 
ground. In both cases the location and extent of damage can 
be relatively accurately detected using the first 5 measured 
natural frequencies and accurately detected using the first 8 
or 10 measured frequencies. 



US 2005/0072234 Al 



17 



Apr. 7, 2005 



[0208] Methods to Handle Some Ill-Conditioned System 

[0209] When the first order perturbation approch is used, 
one needs to solve in each iteration a system of linear 
algebraic equations 

A6G=F (47) 

[0210] which are the linearized equations of (5) and (6), 
where A is the system matrix, F is the vector representing the 
differences between the measured and estimated natural 
frequencies and/or mode shapes, 8G is the optimal changes 
in the stiffness parameters to be found. While the gradient 
and quasi-Newton methods can be used, the generalized 
inverse method is most efficient because it docs not involve 
nested iterations, and thus 8G=A + F, where A + is the gener- 
alized inverse of A. The problem (47) may be ill posed for 
certain cases, especially in the first several iterations, 
because A + can be relatively large and small changes in F 
can result in large changes in the solution 5G. Since the 
stiffness parameters cannot be negative or greater than the 
corresponding values of the undamaged structure, the solu- 
tion may not converge. 

[0211] Ill-conditioning problems do not occur in all the 
examples described above. Sometimes they can occur. Con- 
sider, for example, a cantilever aluminum beam of the same 
dimensions as those of the beam shown in FIG. 12. The 
beam is divided into 36 elements and assumed to have 30% 
of damage at the 5 th element and 90% of damage at the 18 th 
element. The translational degrees of freedom of the first or 
second mode shape vector are used to detect damage. The 
system equations are determinate, and the ill-conditioning 
problem occurs in the iterations. When a natural frequency 
is also used in additional to an incomplete eigenvector 
described above, the ill-conditioning problem can also 
occur. The following regularization methods can be used to 
handle the ill-conditioning problem: 

[0212] Method 1. Estimate A + from (AVl+tiI)- 1 A T , 
where T| is a small positive constant that can be 
searched in several ways. One way is set t|=r|x 
1 .6 lSxmin( 10,||8G|| ,). where , denotes the infin- 
ity norm, with the initial value T|=||A T A|[ 00 , until 
||8G||„ is less than 0.8, for example. 

[0213] Method 2. To constrain the magnitude of 8G, 
we include it in the objective function to be mini- 
mized. For example, we can minimize the following 
objective function 



f(SG) = £ £ (Fi " A « SG ^ + Z {SC>)2 



[0214] instead of (35), where N 0 =n x +n ip N m is the number 
of equations in (5) and (6). The optimal solution can be 
obtained by using the generalized inverse method for the 
expanded system A*=[A;I] and the expanded vector F*= 
[F;0], where I is the mxm identity matrix and 0 is the mxl 
zero vector. 

[0215] Note that the solutions from the regularization 
methods are not strictly the optimal solutions for the original 
objective function in (35). With accurate and sufficient 



measurement information and proper handling of the ill- 
conditioning problem, the system equations can become 
well conditioned in the last few iterations and regulation 
does not need to be applied. Consequently the stiffness 
parameters can be more accurately determined. Sometimes 
regulation may over constrain the magnitude of 8G, and the 
termination criterion ||8G|| co <e may be satisfied during the 
regulation process. In this case, we may set 



)<■ - 3 '' : ' ><;; 



[0216] to amplify the magnitude of 8G, and the iteration is 
terminated after the system equations become well condi- 
tioned. Since the second regulation method does not need to 
search for r\, it can be more effective when it works. When 
one carries out the damage detection procedures using 
different combinations of measured frequencies and mode 
shapes and obtains the same or similar stiffness parameters, 
one can have sufficient confidence on the results obtained. 

[0217] Using only the translational degrees of freedom of 
llic first eigenvector and the first regulation method, the 
estimated bending stiffnesses of all the elements of the beam 
are shown in FIG. 26A. Using the translational degrees of 
freedom of the second eigenvector and the first regulation 
method, the estimated bending stiffnesses of all the elements 
of the beam are shown in FIG. 26B. In both cases the 
locations of damage are exactly detected and the extent is 
relatively accurately determined. Using the translational 
degrees of freedom of the first eigenvector along with the 
first natural frequency and the first regulation method, the 
bending stiffnesses of all the elements of the beam are 
exactly determined, as shown in FIG. 26C. 

[0218] Similarly, using the translational degrees of free- 
dom of the first eigenvector and the second regulation 
method, the estimated bending stiffnesses of all the elements 
of the beam are shown in FIG. 27A. Using the translational 
degrees of freedom of the second eigenvector and the second 
regulation method, the estimated bending stiffnesses of all 
the elements of the beam are shown in FIG. 27B. In both 
cases the locations of damage are exactly detected and the 
extent is relatively accurately determined. Using the trans- 
lational degrees of freedom of the first eigenvector along 
with the first natural frequency and the second regulation 
method, the bending stiffnesses of all the elements of the 
beam are exactly determined, as shown in FIG. 27C. 

[0219] Conclusions 

[0220] Thus, the sensitivities of eigenparameters of all 
orders may, for the first time, be derived using a multiple- 
parameter, general-order perturbation method. The higher- 
order solutions may be used to estimate the changes in the 
eigenparameters with large changes in the stiffness param- 
eters. The perturbation method may be combined with an 
optimization method to form a robust iterative damage 
detection algorithm. The gradient and quasi-Newton meth- 
ods can be used for the first or higher order system equa- 
tions, and the generalized inverse method can be used 
efficiently with the first order system equations because it 
does not involve nested iterations. Including the higher- 
order perturbations can significantly reduce the number of 
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iterations when there is a large level of damage. A modified 
eigenvector expansion method is used to estimate the 
unmeasured component of the measured mode shape. For 
many cases, the location(s) and extent of damage can be 
relatively accurately detected using only measured natural 
frequencies. Methods to handle ill-conditioned system equa- 
tions that may occasionally arise are developed and shown 
to be effective. Numerical simulations on different structures 
including spring-mass systems, beams, lightning masts, and 
frames show that with a small number of measured eigen- 
parameters, the stiffness parameters of the damaged system 
may be accurately identified in all the cases considered. 
Experiments on the different beam test specimens and the 
lightning mast in an electric substation validated the theo- 
retical predictions. The methodology can be readily applied 
to various operation structures of different sizes by incor- 
porating their finite element models or other mathematical 
models. 

[0221] One example of a practical application of this 
damage detection method is shown in FIG. 28, in which a 
structure 10 may represent any type of structure that can be 
modeled by beam elements, on which periodic damage 
assessment must be conducted. This includes, but is not 
limited to, structures such as lightning masts, utility poles, 
cell towers, and other such structures. A structure that can be 
modeled by beam elements is used simply for ease of 
discussion, and it is well understood that the general order 
perturbation method discussed above may be applied to a 
number of different structural members without departing 
from the spirit of the invention as embodied and broadly 
described herein. 

[0222] The structure 10 divided into elements E 1 through 
E n , and with nodes N ± through N n , is equipped with a sensor 
40, such as, for example, an accelerometer. The sensor 40 is 
configured to measure a response to an impact F or a force 
from a shaker applied to the structure 10. The impact F may 
be applied in the form of a single impact, as shown in FIG. 
28, or it may be in the form of a series of random impacts 
F-L through F n , as shown in FIG. 29. The location on the 
beam where the impact F is applied may be varied, but the 
impact F is preferably not applied at one of the nodal points 
of vibration modes whose natural frequencies and mode 
shapes need to be measured. In either instance, the sensor 40 

force from the shaker, and transmits that response to a 
processor 20 equipped with, among other elements, the 
damage detection method 30 discussed above. 

[0223] In accordance with the method as described above, 
upon receipt of the response signal from the sensor 40, the 
processor 20 accesses and performs the method 30 as 
previously discussed, and, as the data converges, determines 
an extent and location, by element, of any structural damage 
present in the structure 10. 

[0224] Although the convergence of the system to reso- 
lution is not dependent on where the impact F is applied to 
the structure 10, a signal to noise ratio of the response may 
be increased, and an efficiency of the system optimized, as 
the application point of the impact F is moved away from a 
fixed point N 0 , if one needs to measure the natural frequen- 
cies and/or mode shapes of lower modes. Similarly, although 
the system will converge regardless of the average energy 
and timing of the impact F or series of impacts F^F^ applied 



to the structure 10, a signal to noise ratio of the response may 
be increased, and an efficiency of the system optimized 
based on a random series of impacts or the random shaker 
excitation to the structure 10, if the structure is large or 
slightly nonlinear or if there is an ambient excitation to the 
structure such as wind. 

[0225] The structure 10 may be excited in a number of 
different manners. For example, an impact F may be applied 
manually or a shaker may be used, at any given point on the 
structure 10, excluding the fixed point N 0 or the nodal points 
of vibration modes whose natural frequencies and/or mode 
shapes need to be measured. However, in an effort to 
increase the signal to noise ratio in the signal provided by the 
sensor 40 to the processor 30 if the structure is large or 
slightly nonlinear or if there is an ambient excitation to the 
structure such as wind and to provide the convenience and 
portability over the shaker test, a series of random impacts 
Fj-F n shown in FIG. 29 may be applied manually or by a 
specially designed dev ice to generate the impact 1 . Although 
accurate and reliable damage assessments can be achieved 
regardless of how the impact F is applied to the structure 10, 
results may be improved using a random series impact 
method, which involves using a series of impacts F of 
random amplitudes and random arrival times. A series of 
random impacts has been shown to increase an energy input 
to the structure 10, improve the signal to noise ratio, 
especially in such situations as strong wind excitation, and 
average out slight nonlinearities that arise, for example, 
from bolted joints and extract linearized eigenparameters. 
[0226] FIG. 29 illustrates a structure 10, sensor 40 and 
processor 20 similar to those shown in FIG. 28. An impact 
F in FIG. 29 is applied to the structure 10, at a location other 
than N 0 or one of the nodal points of vibration modes whose 
natural frequencies and/or mode shapes need to be mea- 
sured, as a series of random impacts F 1 through F n delivered 
at random amplitude and arrival time. These random impacts 
F^F,, may be delivered manually, or, as shown in FIG. 30, 
may he delivered in an automated fashion through the use of 
a random impact hammer device 1550. 
[0227] Different methods have been employed in conven- 
tional vibration testing in order to excite a test specimen. 
Shaker testing, in which a specimen is, simplistically, shaken 
in order to impart a high level of energy, can produce a high 
signal to noise ratio, and can induce random excitation, 
which can average out slight nonlinearities and extract 
linearized eigenparameter parameters. However, shaker test- 
ing is not practically employed in the field on relatively large 
structures, and can be cost prohibitive to conduct. Single 
impact hammer testing addresses the shortfalls of shaker 
testing, in that it is portable and inexpensive to conduct. 
However, single impact hammer testing falls short where 
shaker testing is strong, in that low energy input of single 
impact hammer testing produces a low energy input, a low 
signal to noise ratio with no randomization. To address the 
need for a system which combines the advantages of shaker 
testing and single impact hammer testing, a Random Impact 
Series method for hammer testing is presented which yields 
a high energy, high signal to noise ratio, random system. A 
novel stochastic model is developed to simulate the random 
impact series produced manually and to generate a random 
impact series for a specially designed random impact device. 
[0228] FIG. 31 shows a random impact device 1550, 
according to one embodiment of the invention. Random 
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impact device 1550 includes random signal generating unit 
1560 and a random impact actuator 1570 with impact 
applicator 1580. The impact applicator 1580 preferably 
includes a sensor 1581, such as a force transducer, attached 
at its tip. Sensor 1581 is preferably configured to send data 
to the spectrum analyzer in stiffness parameter unit 103 in 
order to obtain mode shape information, as discussed above 
in connection with FIG. 1A. Sensor 1581 can be coupled to 
the spectrum analyzer in any manner know in the art 
including, but not limited to, wire, optical fiber or even a 
wireless connection. Random impact signal generating unit 
1560 generates outputs 1590 and 1591 to random impact 
actuator 1570. It should be appreciated that outputs 1590 and 
1591 could be coupled to the random impact actuator 1570 
in any manner known in the art including, but not limited, to 
cables, wires, optical fibers, wireless communications, and 
so forth. 

[0229] Output 1590 corresponds to the amplitude as 
discussed below. In particular, random signal generating unit 
1560 outputs a value corresponding to as discussed 
below. Output 1591 corresponds lo T ; , as discussed below. 
[0230] Impact applicator 1580 is shown, for purposes of 
illustration, as impacting structure 1600 in a reciprocating 
motion. Impact applicator 1580 is shown to have an impact 
path 1610, such that when a structure 1600 lies within the 
impact region 1610, impact applicator 1580 impacts struc- 
ture 1600 with a force of random amplitude that arrives at a 
random time, as discussed below. The impact applicator 
1580 and impact region 1610 as shown in FIG. 31 is one 
possible embodiment, and other shapes and impact regions 
may be used while still falling within the scope of the 

[0231] Although random signal generating unit 1560 is 
shown to output two outputs 1590 and 1591 that correspond 
to tjjj and Xj below, it should be understood that signal 
generating unit 1590 could output any signal or signals that 
ultimately results in random impact actuator 1570 driving 
impact applicator 1580 to impact structure 1600 with ran- 
dom arrival times x ; and random amplitudes r|> ; . 

[0232] As discussed above, vibration information is 
obtained by sensor 110, and is sent to stiffness parameter 
unit 103 via sensor coupler 113. The stiffness parameter unit 
103 sends stiffness parameters to the damage information 
processor 123. While the random impact device can be used 
for damage detection purposes as discussed above, it can 
certainly be used just for obtaining the natural frequencies 
and/or mode shapes of the structure for modal testing 
purposes. 

[0233] Experiments conducted on lightning masts by the 
applicant confirm that multiple impact testing performs 
better than single impact testing when there is wind excita- 
tion to the masts. Results are shown for a 65 foot tall mast 
with a 5 foot spike, as shown in FIG. 32, referred to 
hereinafter simply as a seventy foot tall mast. The mast has 
two constant cross section schedule 40 pipes of equal length. 
It can be seen from the results shown in FIGS. 33A-33B that 
the multiple impact test has much better coherence, is less 
noisy away from the resonances, and can pick up the modes 
that were missing in the single impact test. The multiple 
impact test is better at exciting the lower modes that ca n also 
be excited by the wind, which can be seen by comparing the 
frequency response functions (FRF) between 0 and 30 Hz. 



It can been seen from the FRF for the multiple impact tests 
that there are some modes that are missed at near 4 Hz and 
7 Hz with the single impact test, and the mode at 10 Hz is 
much improved compared to the single impact test. 

[0234] A stochastic model will now be discussed which 
describes the random impact series Fj-F^ modeled as a 
Poisson process. The Poisson process is one of a general 
class of processes that arise in problems concerning the 
counting of events in the course of time. The force pulses in 
the series are assumed to have an arbitrary, deterministic 
shape function and random amplitudes and arrival times. 
The force signal in a finite time interval is shown to consist 
of wide sense stationary and non-stationary parts. The 
expectations of the average pow er densities associated with 
the entire force signal and the stationary part of the signal are 
derived and compared, and the power spectral density is 
related to the average power density and the autocorrelation 
function associated with the stationary part of the force 
signal. Numerical simulation is conducted to validate ana- 
lytical predictions, and a relationship between the Fourier 
transform and the discrete Fourier transform used in a 
numerical simulation is produced. Experiments on the four 
bay space frame, as shown in FIG. 10, validated the 
distributions of the random variables and the Poisson pro- 
cess. 

[1(235] Stochastic Model of a Random Impact Series 

[0236] A random impact series is modeled here as a sum 
of force pulses with the same shape and random amplitudes 
and arrival times: 



[0237] where t is time, x(t) is the time function of the force 
signal, x ; e(0, t] is the random arrival time of the i-th pulse, 
and y(t-T;) is the deterministic shape function for all the 
pulses, vj> ; is the random variable describing the amplitude of 
the i-th pulse, and N(t) is the number of the pulses that have 
arrived during the time interval (0,t] and is modeled as a 
Poisson process with stationary increments. All the pulses 
arc assumed to be of width At, and y(t-T,) satisfies y(t-T ; )=0 
if t<t ; and t>x ; +Ax. 

[0238] Since a finite time record is used for the force 
signal in modal testing, we consider only the pulses that 
arrive during the time interval (0,T] of length T, as shown in 
FIG. 34. Note that y(0) and y(Ax) do not have to be equal 
to zero in this model. Because the pulses arriving after time 
t (t<T) will not affect the force signal at time t, the random 
process N(t) in (49) can be replaced by N(T). Also, if a pulse 
arrives at time T, it will vanish at time T+Ax. The last pulse 
in FIG. 34 arrives at time T N(T) and ends at a time between 
T and T+Ax. To include completely this last possible pulse, 
we consider the time interval t e(0, T+Ax]. Equation (49) is 
rewritten as 
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N[T) (50) , . NfT) 

l) = |>*r-T». r.efO.riand re (0,r + Ar] £{| ,,( r -,)}«-. = £ ,.^%( f - ,) e - 



[0239] For the Poisson process N(l) with siaiionary incre- 
ments, the probability of the event {N(t)=n}, where n is an 
integer, is 



[0246] where F denotes Fourier transform, X(jto) is the 
Fourier transform of x(t), to is the angular frequency, and j = 
n^I. Let u=t-T;, then t=u+x ; and dt=du. Equation (54) 
becomes 



[0240] where X is the constant arrival rate of the pulses. By 
replacing t with T, the probability of the event {N(T)=n} is 



[0247] Since the integral in (55) is a deterministic func- 
tion, the expectation of the force spectrum is 



[0241] All the arrival times x i; where i= 1,2,3, . . . , N(T), 
are identically distributed, mutually independent random 
variables. Because the arrival rate of the pulses is constant, 
the arrival limes t ; arc uniformly distributed in [0. T] with 
the probability density function 



[0248] Using the expression for conditional expectations, 
E[E[V|U]]=E(V), where U=N(T) and 



[0242] Similarly, (i-1,2, . . . , N(T)) are identically 
distributed random variables, which are mutually indepen- 
dent and independent (if the distribution of the arrival limes 
X;. While the distribution of r[i ; (i-1,2, . . . , N(T)) is not used 
in the subsequent derivation, it is assumed that i|> ; satisfy the 
Gaussian distribution with the probability density function 



[0249] we have from (56) 

E[XQoj)] = ^'""'] | £ T y(u)e-'"du 



V2?T cr [0250] Since ip ; e (i=l,2, . . . , n) are independent of 
each other and i|\ is independent lift" 1 "" 1 , we have 

[0243] where ^=E[rp ; ] is the mean and ct 2 =E[i|i :1 2 ]-E 2 [i|) ; ] 

is the variance of r[> ; . While the distribution in (53) is not Jy ^.1 = y M = y me -Mn (58) 

used in the analytical derivations, it is used in the numerical [tt ' \~ it it 
simulation and validated experimentally for a manually 
applied random impacts on a four bay space frame, as shown 

in FIG. 10. |- Q251 -| since ^ £ =1 ^ ^ n ) are identically distributed 

[0244] Force Spectrum and its Expectation random variables and so are x ; (i-1,2, . . . , n), we have from 

[0245] The force spectrum is the Fourier transform of the 
force signal in (50): 



XQw) = F[x(t)] = 



(54) 
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[0252] Substituting (59) into (57) yields [0259] Let k=u-v, then u=k+v and du=dk. The double 

integral in (64) becomes 

WWft = (60) 

£ P m (n, T)nE[^jy-^p Tl M^t\£ yWe^du ^ ^ = 

[0253] Substituting (51) and (52) into (60) yields dv£' +&T y(v + k)y(v) e -^ k dk 

E[XQu»] = J[ (Ar) J ] nE ^\l [jy*'" dT ] [0260] Interchanging the order of integration in (65) yields 

J J ),(H)y(v) C -^»-"dHdv = 

£ T y(u)e-*"du J** J*^y(v + *)y(v)«-**«fv 

Jo 

[0261] Let v+k=v, then v=y-k and dy=dv . The first integral 
on the right-hand side of (66) becomes 

[0254] where Taylor expansion of e XT has been used. 
[0255] Average Power Density of x(t) in [0,T+Ax] and its 

Expectation J dk j y(y + k)y(y)e-> uk dv = 

[0256] The average power density of the force signal in r o r ^k 

(50) is defined as J J k ] 0 >^T- k)e ' Md T 



[0262] ('hanging y in (67) back to v and substituting the 

resulting expression into (66) yields 



[0257] where X* (jca) is the complex conjugate of X(jco): 



r jr(u)jr(v)e-W"-» dudv = 



[0258] Substituting (55) and (53) into (62) yields [° 263 1 Notin g tl 



wi r4r (64) r^- 



v(v + H|)v(v)rfv 



[0264] is an even function of k, we have 
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[0265] Substituting (69) into (68) and substituting the 
resulting expression into (64) yields 



> iA m iAiSin[jw(Ti - T m )\ = 0, 0> f (/ = 1,2 N(T)) 



[0268] arc identically distributed random variables, and si 
are x ; (i=l, 2, . . . , N(T)), we have 



E E ^ + E E ^^ cosm ^ - t ») 



l/V(r )l A, : + [/v : (D-/v(r)] 1 A, 2 co- 



[0269] Since X-j and x 2 are independently, uniformly dis- 
tributed random variables in [0,T], the probability density 
function of x 1 -x 2 is 



-continued 



= ATE[i/r 2 ] + 2A 2 £ 2 [^ 
[0271] Substituting (74) into (71) yields 



2A 2 £ 2 [ l A 1 ] 1 " C ° S 2 ( " r) +XTE[^]\ 



[0272] Mean and Autocorrelation Functions of x(t) 

[0273] The first-order cumulant function of x(t) in (50), 
K i[x(t)]> is equal to its mean function, E[x(t)]. The second- 
order cumulant function of x(t), K 2 [x(tj)x(t 2 )], is related to 
its autocorrelation function, E[x(t 1 )x(t 2 )] through 

ElxihXQhKlxitM'^Mk^MQ] (76) 
[0274] where t 1 and t ? are any two time instants in [0,T+ 
Ax]. Following the derivations, the first- and second-order 
cumulant functions of x(t) are 



k- , [ a(/) ]=/■. [ .v(< ) hmViV o T y(t-a)da 
K 2 [^ 1 Wf 2 )]0 Kxx (f 1 , t2 )=XE[ 1 l, 1 2 ]J 0 T y(r 1 -a)y( t2 -a)rfa 



(77) 



[0275] where le[(),T+Ax]. Let l-a=u in (75), then da=-du. 
We have from (77) 



mn\ = \Wx\\ v(«V/« 



[0276] Let 

W(fWoM«)<& (80) 
[0277] Noting that y(u)=0 when u<0 and u>Ax, we have 



from (79) 



[0270] Using (51) and (73) in (72) yields 

Z%i(«.r)(« 2 -«)£[ 1 Af] 



rAE[^]W(0, 0<r<Ar 
E[x(t)] = iAE[^]W(!\T), Lr<t<T 

[\E[l/n][W(&T)-W(t-T)], T < t < T + At 



[0278] where T>Ax is assumed. Let tj-o^u and t 2 -t 1 =k ir 
(78), then da=-du and t 2 -a=u+k. We have from (78) 
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[0279] When |k|>Ax, y(u)y(u+k)=0 and hence k^, t 2 )= 
0. When OSk^Ax, we have from (82) for different t. 



[0284] Fourier Transform of the Mean Function of x( l ) a nd 
its Equivalence to E[X(joj)] 

[0285] Applying the Fourier transform to E[x(t)] in Eq. 
(81) yields 



XE [<p\ ] £ y(u)y(u + k)du, 0 < r, < At - 
£ T k y(u)y(u + k)du, AT-k< tl < 
XEW\]^\{u)y(u + k)du, T < h < T + . 



[0280] When -Ax£k<0, we have from (82) for different t 2 

, fi (84) 
XEWl ] J y(u)y(u + k)du, -tsi,< Ar 

AEWr? ] J^MM" + k)du, Ar< h <T-k 

XE[^]J^y(u)y(u + k)du, T-k< h <T + Ar 



[0281] Let u+k=v in (84), then u=v-k and du=dv. We have 
from (84) after changing v back to u 



XEi^P^ y{u-k)y(u)du, 
A£[^]j^ A " y(u-k)y(u)du, 
XE[^]j^y(, 



y(u-k)y(u)du, T — 



[0282] Combining the second equations in (83) and (85), 
we have for t t and t 2 in [Ax,T] 



[0283] Since by the second equation in Eq. (81), E[x(t)] is 
a constant for t e[Ax, T], and by (86), K xx (t 17 tj is a function 
of k=t 2 -t 1 for t-L and t 2 in [Ax, T], x(t) is a wide-sense 
stationary random process in [Ax, T]. Substituting the sec- 
ond equation in (81) and (86) into (76) yields the autocor- 
relation function for t t and t 2 in [Ax,T] 



r 4 v 



|ll (Ar)-M (t-T)-\e->"dt\ 



[0286] The three integrals in (88) are referred to as I 1; I 2 , 
and I 3 , respectively. Consider first the third integral in (88) 



W?l J*" m yWy(ui 



[0287] Let t-T=6 in (88), then t=T+6 and dt=d6. We have 
from Eq. (89) 

7 3 =e-i^ 0 ^W(AT)e- i " e de-e- i ° )T; 0 A, W(e)e- i ° )e de (90) 
[0288] Changing 0 in (90) back to t and combining I , and 
I 3 yields 

/ 1 +/ 3 =e- i ° )T ^ 0 A, W(AT)e- i ° ,, A+(l-e- i ° , ' I )J 0 Ai; ty(t)e- 

j, *dt (91) 

[0289] Adding I, to (91), simplifying the expression, and 
substituting it into (88) yields 



F \E[x(t)]} = ( 
e~'"dr + j tV(Ar)e->"*£fr] = A£W I ] 



[0290] We will show that F{E[x(t)]j in (92) is equivalent 
to E[X(jw)] in (91). By (80), we have 



I(t)dt, re[0,Ar] 



[0291] The integral in (61) can be written as 

( 

-£ T y M du-J.£ T y M £e-^d V du 
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[0292] Interchanging the order of integration in the double [0301] where X s (j(fl)=| -Vr r x(t)e JOJt dt. Taking the expec- 
integral in (94) and using (93) yields tation of (96) yields 

£ T y W e-^d U = <95) EfrW-T^fc^EW^V-**!^** ^ 

£ y(u)du-j«>f o dv£ y(u)e-^du ^ k=t 2 -t t in (97), then dk=dt 2 . Since 

E[x(t 1 )x(t,)]=R^(k), (97) becomes after interchanging the 
= W(Ar)|l + jojj' T [^- - l] c -^dv} order of integration 



[0293] Substituting (95) into (61) yields (91). This shows 
that the expectation E and the Fourier transform F are 
commutative as both are linear operators. 

[0294] Equation (92) consists of two parts: the first part, 



AE[lAl](l-f~""')W(Ar)^- 



=0-4-^)-" 

[0303] Substituting (87) into (98) yields 
E[S 2 M]=AEWfi 



[0295] is the Fourier transform of the stationary part of 
E[x(t)] in [Ax,T], and the second part, 2A 2 £ 2 [iA,]w 2 (Ar)J^ ~ Ar - J^LjeMj k 

AE[(W(i -e '""lU'tAri - i)« "-"dr. [0304] Noting that y(u)y(u+|k|)=0 when |k|>Ax and 

[0296] is the sum of the Fourier transforms of the nonsta- X yb + W>yW>^ - 

tionary parts of E[x(t)] in [0,Ax] and [T,T+Ax]. When Ax^O, 

[0305] an even function of |k|, we have from (99) 



since 



e[5;( (U )] = :a : /; : [ia,|u- : ( a 



[0297] is finite, 



"')'(« + \k\)y{u)du[\ - 



[T - At) 

icosukdk 



[0306] where T>2Ax has been assumed. 

[0307] The power spectral density of x(t) can be obtained 

from (100) by increasing T to infinity 



[0298] and consequently the second part of E[X(jco))], 

approaches zero. 

[0299] Average Power Density of x(t) in [Ax,T], Its Expec- A£[l/, ' ] £K ' ' ( ' 

tation, and Power Spectral Density 

[0300] Since x(t) is stationary in [Ax,T], the average power [ 0308 ] where we have used 
density of x(t) in [Ax,T] is defined as 



5 2 (w) = 2^A 2 £ 2 [iAi]M /2 (At) i 5M+ (101) 
\k\)I(u)cos(u>k)dk 
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[0309] in which 8(0) is the Dirac delta function. The 
power spectral density can also be obtained from (98) by 
increasing T to infinity 



.S>„=li m T— l-.|S.,„„|= 
.■RJkv '■\lk=l']RJk>\ 



(103) 



[0310] where R xx (-k)=R xx (k) has been used. Equation 
( 1 03) is (he well-known Wiener-Khintchine theorem, which 
states the power spectral density is the Fourier transform of 
the autocorrelation function. Substituting (87) into (83), and 
noting that I(u)I(u+|k|)=0 when |k|>Ax and the Fourier 
transform of 1 is 2;to((>>), yields (101). Note that the power 
spectral density is only defined for a wide -sense stationary 
process with an infinite time record. When the mean ampli- 
tude of each pulse EfipJ is not equal to zero, there is an 
associated delta function in the power spectral density. 

[0311] Comparison of EfS^w))] and E[S 2 (w)] 

[0312] By (100), E[S 2 (w)] consists of two parts: the first 
part, 



[0317] When y(v)=W 0 8(v), where W 0 
have from (104) 



v!v + |A-|lY!vir/.'cos„,/u/A = 



Hm|j^ A V 0( 5(v)£-^dv| =W 0 2 



[0318] Since |k|<Ax in (100), we have |k|-»0 as Ax^>0 and 
hence 



[0313] which depends on the arrival rate X, the mean 
amplitude of each pulse E^-J, the total area of the normal- 
ized shape function W(Ax), and the time length T-Ax, 
describes the average effects of the stationary part of x(t) in 
[Ax,T] and is referred to here as the first-order statistical 
power density, and the second part, 



[0319] Substituting ( 1 05) and ( 1 06) into ( 1 00) and noting 
that W(Ax)=W 0 when y(v)=W 0 8(v), and substituting (105) 
into (75) and noting that T+Ax-»T as Ax-»0, yields 



AE[iA?]J^ AT m y(u + \k\)y(u)di\l - JL_jco«a*<f, 



£[5 1 M] = 2A 2 £ 2 [iA I ]- 



[0314] which depends on the mean square amplitude of 



ary part of x(t) and is referred to here as the second-order 
statistical power density. While E[S 1 (co)] in (75) consists of 
two parts, the first part, 



v,v + H|,v,v„/,vo«,„A «lk\- 



[0315] depends also on the shape function y(0) as the 
shape function is used in calculating the power associated 
with the nonstationary parts of x(t) in [0,Ax] and [T, T+Ax]. 

[0316] When the shape function is a delta function (i.e., 
Ax-»0), the nonstationary parts of x(t) vanish and E[S,((n)] 
in (75) can be shown to be equivalent to E[S,(<.>)] in (100). 
By (68) and (69), we have 



[0320] By (104) the power spectral density in (101) can 1 
expressed as 

S»=2 3t X, 2 £ 2 [^ 1 ]W 2 (AT)8((o)+A£[ 1 l. 1 2 ]|rO-»)| 2 (10 

[0321] where Y(j w)=J 0 AT y(t)e- jrot . When T^oo in (75), v 
have by using (104) 



[0322] Equation (109) has a slightly different form from 
that of (108) because the power associated with the nonsta- 
tionary parts of x(t) is included in (109). When y(t)=W 0 8(t) 
(108) and (109) reduce to 

5»=S((0)=2 3t ?>, 2 ir 2 [^ 1 ]W 0 2 5((0) + A£[^ 1 2 ]W 0 2 (110) 

[0323] Equation (110) can also be obtained from (106) by 
letting T-»o° and using 
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EXAMPLES AND NUMERICAL SIMULATION 

[0324] When the shape function of the pulses is repre- 
sented by a half sine wave, i.e., 



[0325] where H(0) is the Heaviside function, we obtain by 
using (60), (64), and (103) 



shown as a solid line in FIG. 38, is calculated using (75) . 
The curve for 10 logjEfS^ffi)]} with X=l/s and the other 
parameters unchanged is shown as a dashed line in FIG. 38. 
It is seen that Ep^co)] increases by 4.14 to 15.6247 times 
in the frequency range shown when X is increased from 
X^l/s to X 1 =4.14/s. This result can be shown by using (75) 



lim£[5,(^)]| i=l2 



Ai/.-|,.v,|/- + 



//-hi = 
■://.i-ni 



E[X(p,)] = 



-2 A - /_-.<> i »o- ./-I. I 



E[S 2 (co)]=XE(^)At 2 - 



(w 2 Ar 2 -7T 2 ) 3 (7--Ai-) 



[+(8sin(wAr)w^ + 
27cos(o;At>V +27>V)Ar 2 ] 



[0327] This shows that a larger arrival rate X would 
increase the energy input to the structure over the entire 
frequency domain. 

[0328] Numerical simulation is undertaken next to vali- 
date the analytical predictions. The random number N(T) 
satisfying the Poisson distribution in (5 I ) w ith >.=4. 14 s and 
T=8 s is generated using MAT! AH. Similarly, the random 
numbers corresponding to the random variables %■ (i=l, 2, . 
. . , N(T)), satisfying the uniform distribution in (52), and the 
random numbers corresponding to the random variables t|) ; 
(i=l, 2, . . . , N(T)), satisfying the Gaussian distribution in 
(53) with ,m=0.8239 N and a 2 =0.7163-0.82392 N 2 =0.0375 
N 2 , are generated. Using the shape function constructed 
earlier, a sample function of x(t) in (50) at time t=rh, where 
r=0, 1, , 



8AA." : i w i ,Ar : (! - co<**T - 



[0326] Consider next the normalized shape function y(t) 
shown in FIG. 35 with unit maximum amplitude. It is 
obtained by averaging a series of normalized force pulses 
from impact tests on the four-bay space frame as shown in 
FIG. 10. There are 21 sample points in the shape function, 
which are connected, as shown in FIG. 34. Other parameters 
used are T=8 s, Ax=20xT/l 024=0. 15625 s where h=T/1024= 
0.0078 125/s is the sampling interval, X=4.14/s, Ef^^ 
0.8239N, and E[i|- , 2 ]=0.7 1 63N 2 . The curve lor E[x(t)] in the 
time interval from 0 to 8.15625 s, shown as a solid line in 
FIG. 36, is calculated using (81). It is seen that E[x(t)] 
increases from 0 to 0.0352 N in the first 0.15625 s, remains 
at 0.0352 N from 0.15625 s to 8 s, and decreases from 
0.0352 N to 0 in the last 0.15625 s. The curve for 20 
log|E[X(jo>)] in the frequency range from 0 to 50 Hz, shown 
as a dotted line in FIG. 37, is calculated using (61). The 
curve for 10 loglEfSjQff)))]} in the same frequency range, 



[0329] denoted by x r , can be obtained. The discrete Fou- 
rier transform (DFT) of the time series {x r } is calculated 
using MATLAB. The DFT of the series {x r } is defined by 



-it"* 



[0330] where q=0, 1, . . . , R-l. Equation (118) is an 
approximate formula for calculating the coefficients of the 
Fourier series of a periodic function whose values in the 
period [0, T+Ax] are given by those of x(t): 
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[0331] The Fourier components X q correspond to harmon- 
ics of frequency 



[0332] Recall that the Fourier transform of x(t) in (50) is 

given by 



X(Jo) = £ +&T x(t)e->" J < dt 



[0333] Let 



W "« T + At 



[0334] in (120) and compare the resulting expression with 
(119), we find that X q in (118) multiplied by T+Ax provides 
an approximate value of X(jco) at frequency co q . Similarly, 
X X* q multiplied by T+Ax provides an approximate value of 



X(ju>)X'(ju>) 



[0335] at frequency co q . By averaging 5000 sample series 
of {x r }, we obtain the curve for E[x(t)], shown as a dotted 
line in FIG. 36. By averaging 1000 sample series of {20 
log[(T+AT)|Xj]). we obtain Hie curve for 20 log|L[\(j,.,)]|. 
shown as a solid line in FIG. 37. By averaging 100 sample 
series of {10 log[(T+Ax)|X q X* q |]}, we obtain the curve for 
10 logjE^Qo))]}, shown as a dotted line in FIG. 38. The 
numerical results are in good agreement with the analytical 

[0336] The stochastic model was experimentally validated 
for an experimenter conducting manually a random series of 
impacts on the four bay space frame as shown in FIG. 10. 
One hundred ensemble averages were used. The experimen- 
tal probability density functions of the arrival time, the 
number of arrived pulses, and the pulse amplitudes are in 
good agreement with the analytical values, as shown in 
FIGS. 39-41. 

[0337] Thus, the system and method for detecting struc- 
tural damage and the random impact series method as 
embodied and broadly described herein can be applied to an 
unlimited number and type of structures to provide auto- 
mated, reliable damage detection and assessment and to 
conduct modal testing. This system could be further auto- 
mated to conduct periodic tests and provide results to a 
centralized monitoring section. Regular health monitoring of 
these types of structures could provide additional protection 
against potential failure, as well as a characterization of 
usage and wear over time in particular environmental con- 
ditions for predicting useful service life. 
[0338] The foregoing embodiments and advantages are 
merely exemplary and are not to be construed as limiting the 



present invention. The present teaching can be readily 
applied to other types of systems. The description of the 
present invention is intended to be illustrative, and not to 
limit the scope of the claims. Many alternatives, modifica- 
tions, and variations will be apparent to those skilled in the 
art. In the claims, means-plus-function clauses are intended 
to cover the structures described herein as performing the 
recited function and not only structural equivalents but also 
equivalent structures. 
What is claimed is: 

1. A system for determining stiffness parameters of a 
structure, comprising: 

a sensor arranged to measure vibrations of said structure 
and output vibration information; and 

a stiffness parameter unit for receiving said vibration 
information, determining natural frequency data of said 
structure, and determining the stiffness parameters of 
said structure using said natural frequency data. 

2. The system according to claim 1, further comprising 
multiple sensors arranged to measure vibrations of said 
structure and output vibration information. 

3. The system according to claim 1, wherein said stiffness 
parameter unit comprises an iterative processing unit. 

4. The system according to claim 1, wherein said stiffness 
parameter unit comprises an outer iterative processing unit 
and an inner iterative processing unit. 

5. The system according to claim 3, wherein said iterative 
processing unit determines said stiffness parameters using a 
first order perturbation process. 

6. The system according to claim 3, wherein said iterative 
processing unit determines said stiffness parameters using a 
higher order perturbation process. 

7. A system for determining stiffness parameters of a 
structure, comprising: 

a sensor arranged to measure vibrations of said structure 

and output vibration information; and 
a stiffness parameter unit for receiving said vibration 

information and determining said stiffness parameters 

with an iterative processing unit. 

8. The system according to claim 7, wherein said iterative 
processing unit comprises an outer iterative processing unit 
and an inner iterative processing unit. 

9. The system according to claim 7, wherein said iterative 
processing unit determines said stiffness parameters using a 
first order perturbation process. 

10. The system according to claim 7, wherein said itera- 
tive processing unit determines said stiffness parameters 
using a higher order perturbation process. 

11. A stiffness parameter unit for determining stiffness 
parameters for a structure, comprising: 

an input for receiving vibration data related to the struc- 

an analyzer for converting said vibration data to spectral 
data; and 

an interative processing unit for receiving said spectral 
data and outputting said stiffness parameters using 
natural frequencies of the structure. 

12. A stiffness parameter unit for determining stiffness 
parameters for a structure, comprising: 

an input for receiving vibration data related to the struc- 
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an analyzer for converting said vibration data to spectral 
data; and 

an interative processing unit for receiving said spectral 
data and outputting said stiffness parameters using a 
perturbation process. 

13. The stiffness parameter unit according to claim 12, 
wherein said perturbation process comprises a first order 
perturbation process. 

14. The stiffness parameter unit according to claim 12, 
wherein said perturbation process comprises a higher order 
perturbation process. 

15. A system for determining damage information of a 
structure, comprising: 

a sensor arranged to measure vibrations of said structure 
and output vibration information; 

a stiffness parameter unit for receiving said vibration 
information, determining natural frequency data of said 
structure, and determining the stiffness parameters of 
said structure using said natural frequency data; and 

a damage information processor for receiving said stiff- 
ness parameters and outputting damage information. 

16. The system according to claim 15, wherein said 
damage information processor outputs damage location 
information or extent of damage information. 

17. A system, comprising: 
a structure; 

a sensor arranged to measure vibrations of said structure 
and output vibration information; and 

a stiffness parameter unit for receiving said vibration 
information, determining natural frequency data of said 
structure, and determining the stiffness parameters of 
said structure using said natural frequency data. 

18. The system according to claim 17, further comprising 
a damage information processor for receiving said stiffness 
parameters and outputting location of damage. 

19. The system according to claim 18, wherein said 
damage information processor comprises a damage location 
processor for determining damage location information. 

20. The system according to claim 18, wherein said 
damage information processor comprises a damage extent 
processor for determining extent of damage information. 

21. The system according to claim 18, wherein said 
damage information processor comprises a damage extent 
processor for determining extent of damage information and 
a damage location processor for determining damage loca- 
tion information. 

22. The system according to claim 17, wherein said sensor 
comprises a velocimeter. 

23. The system according to claim 17, wherein said sensor 
is attached to said structure. 

24. The system according to claim 17, wherein said sensor 
is not attached to said structure. 

25. The system according to claim 17, wherein said 
stiffness parameter unit further comprises a spectral ana- 
lyzer. 

26. The system according to claim 17, wherein said 
structure comprises a beam. 

27. The system according to claim 17, wherein said 
structure comprises a truss. 

28. The system according to claim 17, wherein said 
structure has a longest dimension less than 1.5 meters. 



29. The system according to claim 17, wherein said 
structure has a longest dimension less than 2.5 meters. 

30. The system according to claim 17, wherein said 
structure has a longest dimension less than 10 meters. 

31. The system according to claim 17, wherein said 
structure has a longest dimension less than 50 meters. 

32. A device, comprising: 

a random signal generating unit for generating first and 
second outputs; 

a random impact actuator for receiving said first and 
second outputs; and 

an impact applicator coupled to said random impact 
actuator and having an impact region; 

wherein said random impact actuator drives said impact 
applicator such that the force and arrival times of said 
impact applicator at said impact region are random. 

33. The device of claim 32, wherein said random impact 
actuator drives said impact applicator in accordance with 
said first and second outputs. 

34. The device of claim 33, wherein the first and second 
outputs comprise independent random variables. 

35. The device of claim 34, wherein the first and second 
outputs determine the force and arrival times, respectively, 
of the impact applicator at the impact region. 

36. A system, comprising: 

a structure; 

a random impact device for inducing vibrations in said 
structure; 

a sensor arranged to measure vibrations of said structure 
and output vibration information; and 

a stiffness parameter unit for receiving said vibration 
information, determining natural frequency data of said 
structure, and determining the stiffness parameters of 
said structure using said natural frequency data. 

37. The system of claim 36, wherein the random impact 
device comprises: 

a random signal generating unit for generating first and 
second outputs; 

a random impact actuator for receiving said first and 
second outputs; and 

an impact applicator coupled to said random impact 
actuator and having an impact region; 

wherein said random impact actuator drives said impact 
applicator such that the force and arrival times of said 
impact applicator at said impact region are random. 

38. The device of claim 37, wherein said random impact 
actuator drives said impact applicator in accordance with 
said first and second outputs. 

39. The device of claim 38, wherein the first and second 
outputs comprise independent random variables. 

40. The device of claim 39, wherein the first and second 
outputs determine the force and arrival times, respectively, 
of the impact applicator at the impact region. 

41. A system for determining stiffness parameters of a 
structure, comprising: 
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a sensor arranged to measure vibrations of said structure 
and output vibration information; and 

a stiffness parameter unit for receiving said vibration 
information, determining mode shape information, and 
determining the stiffness parameters of said structure 
using said mode shape information. 

42. The system according to claim 41, further comprising 
multiple sensors arranged to measure vibrations of said 
structure and output vibration information. 

43. The system according to claim 41, wherein said 
stiffness parameter unit comprises an iterative processing 
unit. 



44. The system according to claim 41, wherein said 
stiffness parameter unit comprises an outer iterative process- 
ing unit and an inner iterative processing unit. 

45. The system according to claim 43, wherein said 
iterative processing unit determines said stiffness parameters 
using a first order perturbation process. 

46. The system according to claim 43, wherein said 
iterative processing unit determines said stiffness parameters 
using a higher order perturbation process. 



